magattacaのブログ

日付以外誤報

OpenMMをステップバイステップで 〜Part 6 MDAnalysisでトラジェクトリ解析(①)〜

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

前回まででOpenMMを使ったMDシミュレーション準備&実行(とコードの書き方の理解)ができました。今回はシミュレーション結果のトラジェクトリ解析に挑戦したいと思います。

この記事で取り上げるソフトウェアはMDAnalysisです。またコードはTeach OpenCADD T20・Analyzing molecular dynamics simulationsを参考にしています。

途中で力尽きたので今回は解析のパート①です。*1

MDAnalysisでゆらゆら解析!

1. MD計算のデータ解析って何するの?

そもそもMD計算のデータ解析はどんなことをすれば良いのでしょうか?

私のような初学者にとって素晴らしい講演動画と資料がありました!

高度情報科学技術研究機構(RIST)のオンラインサロン 第4回 スパコンコロキウム でおこなわれた埼玉大学 松永 康佑先生のご講演です(研究室HP)。


www.youtube.com

発表資料(PDF):生体分子系の分子動力学シミュレーションデータの解析入門

基本的な考え方から、②データ解析入門、③携わられたご研究の実例まで紹介してくださっていてとても楽しい講義です。

特に印象に残った点は「シミュレーション」(データの生成)と「データ解析」の「分業」に分野の研究スタイルが変化しているということです。専用のスパコンで計算したデータが公開され、別のグループがそのデータを使って新規な解析手法の論文化を行う例が増えているそう。オープンサイエンスの力を感じますね。

また、スライドp19「データに溺れないために:データ解析の鉄則」はMDシミュレーションの解析に限らない大事な考え方のように思いました。

紹介されていた研究の実例はこちら。薬剤耐性化に関与する多剤排出トランスポーター AcrBの機構の解明です。

薬剤排出メカニズムのMDシミュレーション解析例

  • 論文:Y. Matsunaga, T. Yamane, T. Terada, K. Moritsugu, H. Fujisaki, S. Murakami, M. Ikeguchi, and A. Kidera,
    “Energetics and conformational pathways of functional rotation in the multidrug transporter AcrB”,
    eLife 7, e31715 (19 pages) (2018) DOI: 10.7554/eLife.31715
  • 日本語解説:「全原子分子動力学シミュレーションが解き明かす多剤排出トランスポーターAcrBの薬剤排出機構」 生物物理 2019(59)084-087

(小さな)プロトン化状態の変化で(大きな)膜貫通部位の運動がおこるというのは驚きです。タンパク質は本当によくできていますね。

2. この記事で試す解析とパッケージ

タンパク質の動的な挙動の理解はとても魅力的ですが、とても素人に手を出せるものではないので、ここでは基本的な解析手法の動かし方を試すことにします。

具体的には以前このブログでも取り上げたMaking-it-rainを参考に以下の解析を行います。

また利用するMDトラジェクトリ解析パッケージはMDAnalysisです。

MDAnalysisはTeachOpenCADDT20・Analyzing molecular dynamics simulationsでも使われていて、Pythonから利用できます。また、チュートリアルが充実かつ開発コミュニティが活発そうな感じなのも魅力的です。あとロゴが可愛い。

インストールはpipcondaでOKです(MDAnalysis - Installation)。

mamba install mdanalysis -c conda-forge 

また、トラジェクトリの可視化を行うためにNGLviewを利用します(GitHub, Document)。

こちらもpipcondaでインストールできます。

mamba install nglview -c conda-forge

mamba速い。

3. シミュレーション全体の様子を眺めよう

3-1. NGLviewerで動画を見よう

冒頭に挙げた松永先生のコロキウムでは、まず最初に「シミュレーションがうまくいったか?動画で眺めること。」を推奨されていました。失敗しているとタンパク質が崩壊するそうです。

動画で確認するために、MDAnalysisを利用してシミュレーション結果のファイルを読み込みましょう。前回の記事(Part 4)で得られたファイルは、「① OpenMMReporterを使って出力したトラジェクトリのDCDファイル」と「②PDBファイル形式の前処理後の系の構造」です。

MDAnalysisはUniverseと呼ばれるオブジェクトを中心のデータ構造としてさまざまな処理を行う構成になっているようです(MDAnalysis tutorial - Basics)。

Universeオブジェクトはtopologytrajectoryを読み込んで構築できます。

import MDAnalysis as mda

# topologyとtrajectoryからUniverseオブジェクトを構築
md_universe = mda.Universe("simulation_data/3poz-processed.pdb", "simulation_data/trajectory.dcd")

print(type(md_universe))
#  <class 'MDAnalysis.core.universe.Universe'>

次に、Jupyter notebook上で動画を確認するためのツールとしてNGL viewernglviewを利用します。

nglviewにはさまざまなファイルを可視化する機能があらかじめ用意されており、show_ファイル形式みたいな感じで使えます(nglview : Usage)。

MDAnalysysを使うならshow_mdanalysys()として、Universeオブジェクトを渡せばOKです。

import nglview as nv

view = nv.show_mdanalysis(md_universe)
view

ぱっと見、大きな崩壊はなさそうですね!よかった。

自動で水やイオンを消して見やすくしてくれるので便利です。

3-2. NGLviewが表示されない?

もしviewを実行してもjupyter notebookに何も表示されない場合は以下をお試しください。

jupyter notebookのバージョンがv5.0の時は原子数の多い大きな系入力できなくなっている可能性があるので、ノートブック立ち上げ時にフラグを追加して上限を増やしてあげます。

jupyter notebook --NotebookApp.iopub_data_rate_limit=10000000

もしくはwidgetnbextensionnglviewをノートブックで使えるようにアクティベートする必要があるかもしれません。

NGLviewのGitHub READMEには以下のコマンドが書かれています。*2 *3*4

python -m ipykernel install --sys-prefix
jupyter nbextension enable --py --sys-prefix widgetsnbextension
jupyter nbextension enable --py --sys-prefix nglview

私の環境では以上を実行すると表示されるようになりました。

3-3. エネルギー変化のプロット

今回のシミュレーションではポテンシャルエネルギー(kJ/mol)や温度(K)の継時的な変化を「log.txt」ファイルに書き出しました。

具体的な解析の前にこちらも確認しておきましょう。タブ区切りなのでPandasのread_table()でOKです。

import pandas as pd

log_df = pd.read_table("simulation_data/log.txt")
log_df.head()

import matplotlib.pyplot as plt

x = log_df["Step"]
y = log_df["Potential Energy (kJ/mole)"]

plt.plot(x, y)

plt.title("change of the potential energy")
plt.xlabel("steps (per 0.002ps)")
plt.ylabel("potential energy [kJ/mol]")
plt.show()

横軸がStepで縦軸が系のポテンシャルエネルギーです。一定の範囲内で揺れ動いているようです。

今回のシミュレーションではStep Sizeを「0.002 ps」として、「500,000 steps」計算しているので1 ナノ秒のシミュレーションです。1000 steps毎にログを出力しているのでデータポイント2 ps毎の500個になります。

温度が保たれているかプロットしておきましょう。

x = log_df["Step"]
y2 = log_df["Temperature (K)"]

plt.plot(x, y2)

plt.title("change of the temperature")
plt.xlabel("steps (per 0.002ps)")
plt.ylabel("temperature [K]")
plt.show()

設定の300Kから-4K ~ +4K程度の範囲内に収まっていました。これは十分に温度一定のアンサンブルとみなして良いのでしょうか??

よくわからないので有識者の方ご教示いただけると嬉しいです。

4. MDAnalysisでトラジェクトリ解析

なんとなく上手くいってそう(平衡状態のサンプリングができていそう)なので、MDAnalysisを使ったより具体的な解析に進みます!

4-1. MDAnalysisのオブジェクトと解析対象の選択

まず最初にMDAnaysisにおいて一番基本的なオブジェクトについて触れておきます。

先に述べた通りMDAnalysisでは「topologytrajectoryから構築するUniverseオブジェクト」に対してさまざまな操作を行うことで解析を行います。

Universeの重要な属性の一つがUniverse.atomsです。これはMDAnalysisの「最も基本的な構成要素であるAtomオブジェクト」のリスト(AtomGroup)です。読んで字のごとく系を構成する原子ですね。

例えば「①Universeオブジェクト、② atomsの最初の一つ、③ atomsの最後の一つ」を取り出してみるとこんな感じです。

print("whole system : ", md_universe)
print("example atoms : ", md_universe.atoms[0])
print("example atoms : ", md_universe.atoms[-1])

# whole system :  <Universe with 57648 atoms>
# example atoms :  <Atom 1: N of type N of resname GLN, resid 701 and segid A and altLoc >
# example atoms :  <Atom 57648: Cl of type Cl of resname CL, resid 20103 and segid 4 and altLoc >

「① 系全部で57648原子」あり、その中には「②タンパク質の701番残基GLN窒素原子(N)」や「③加えたイオンの塩素原子(CL) 」があることがわかります。

MDシミュレーションのトラジェクトリ解析では「のうち特定の要素(タンパク質の特定の残基、原子、リガンド)に焦点を当てて時間変化を追跡したい!」といったことが多いようです。

このような時にはMDAnalysis.Universe.select_atoms()メソッドを使うことで、引数で指定した条件のAtomGroupを選択して解析を行うことができます。

例えば「① タンパク質 ("protein")、② 主鎖("backbone")、③ Cα原子 ("name CA")」を選択したい場合は以下のようになります。

print("protein atoms : ", len(md_universe.select_atoms("protein")))
print("backbone atoms : ", len(md_universe.select_atoms("backbone")))
print("CA atoms : ", len(md_universe.select_atoms("name CA")))

# protein atoms :  5120
# backbone atoms :  1268
# CA atoms :  317

それぞれ選択した後lenatomの数を数えています。「②主鎖の原子数」が「③Cα原子」のちょうど4倍になっているので正しく選択できていそうです。

今回対象としているタンパク質(EGFR)のデータは「PDB ID : 3POZ」です。上図の通り前処理の段階でN末端とC末端を削っているので、シミュレーション構造には317残基残っています。Cα原子の数に一致しているようで良かった!

select_atoms()の引数には他にも色々な選択肢があるようです。ドキュメントの3. Selection commandsなどをご参照ください。

4-2. アラインメント

MDAnalysisの基本的な構成がわかったところで、解析の準備段階として、トラジェクトリにおけるタンパク質のアラインメントをとります。

TeacOpenCADD T20によると、(今回のように)シミュレーション時間が短い時はあまりおこらないそうですが、タンパク質が回転したり移動したりすることがあるそうです。大きく位置がずれると解析が難しくなるのであらかじめアラインメントをとって揃えておく方が良いそうです。

また、MDAnalysisの「ドキュメント: Alignments and RMSfitting」にも、「RMSDRMSFといった値が意味のあるものにするためにはトラジェクトリのアラインメンとを取ること」との記述があります。

アラインメントにはMDAnalysis.align.AlignTrajを使います(Aligning a trajectory to a reference)。ここではトラジェクトリの最初の構造リファレンスとして全体を合わせます。

from MDAnalysis.analysis import align

# トラジェクトリの位置を最初のフレームにリセットする
md_universe.trajectory[0]

# align.AlignTrajでアラインメントを取る
alignment = align.AlignTraj( mobile = md_universe, reference = md_universe, select="protein", in_memory=True)
alignment.run()

動画で見ます。

view = nv.show_mdanalysis(md_universe)
view

・・・違いがわからない。

AlingTrajの引数ではmobileで動かす対象のUniverseオブジェクトを指定し、referenceUniverseオブジェクトに合わせます。上記ではどちらも同じUniverseオブジェクトを指定してますが、実行前にトラジェクトリの最初のフレームにリセットする操作を行なっているため、開始時点の位置に対して全体のアラインメントをとることになっています。

また、引数selectではトラジェクトリのうち「どの構造のアラインメントを取るか」選択することができます。上記では「"protein"」とすることでタンパク質全体を選んでいます。先に見たように「"name CA"」でCA原子(アミノ酸Cα炭素原子)を選択することもできます。

最後に、今回は引数「in_memory=True」とすることでその場でアラインメントをとったものに置き換えていますが、トラジェクトリのサイズが大きくメモリに載せきれない場合は、一度ファイルに書き出し再度トラジェクトリを読み込むといった処理もできます。

以下のようにすればOKです。

# 検証用にuniverseオブジェクトを作成
md_universe2 = mda.Universe("simulation_data/3poz-processed.pdb", "simulation_data/trajectory.dcd")
md_universe2.trajectory[0]  

# アラインメントをファイルに書き出す場合。
alignment2 = align.AlignTraj(mobile = md_universe2,  reference = md_universe2,  select='name CA', filename='simulation_data/CA_aligned.dcd',  in_memory=False)
alignment2.run()

CA原子でアラインメントをとって「CA_alinged.dcd」ファイルに書き出しました。*5読み込んで可視化します。

md_universe3 = mda.Universe("simulation_data/3poz-processed.pdb", "simulation_data/CA_aligned.dcd")

view = nv.show_mdanalysis(md_universe3)
view

動画のどこが違うがお分かりでしょうか?・・・私にはさっぱりわかりません。

4-3. RMSD

アラインメントも無事とれたので、とっつきやすいRMSDwikipedia - 原子位置の平均二乗偏差)から解析を行なってみましょう!

だいたいこういう話です。*6

少し注意点ですが、MDAnalysisのメソッドには小文字のrmsd大文字のRMSDがあり、解析機能(や引数)が異なります。

4-3-1. フレーム2点間のRMSD(rms.rmsd()

まず、RMSDの評価としてトラジェクトリの「フレームの最初最後でどれだけ構造変化があったか?」を検証してみます。2つの状態の比較の場合はMDAnalysis.analysis.rms.rmsd()を使えば良いようです。rmsd小文字であることにご注意ください(MDAnalysis Tutorial : 8. Using the MDAnalysis.analysis modules)。

from MDAnalysis.analysis import rms

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

# 最初のタンパク質の位置
protein_first = md_universe.select_atoms("protein").positions

# 最後のフレームに移動
md_universe.trajectory[-1]

# 最後のタンパク質の位置
protein_last = md_universe.select_atoms("protein").positions

# rmsdを評価
print(rms.rmsd(protein_first, protein_last))

#   1.7643614273332058

1.76Å・・・なるほど?・・・この値がどの程度なのか良くわかりません。

4-3-2. トラジェクトリ全体のRMSD推移(rms.RMSD()

次にタンパク質の主鎖backbone)に焦点を当てて、トラジェクトリ全体での変化を見てみましょう。全体を見る場合はMDAnalysis.analysis.rms.RMSD()を使います(大文字RMSDに注意!)。

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

# RMSD解析
RMSD_analysis = rms.RMSD(md_universe,   # アラインメントを取るオブジェクト  
                     md_universe,     # アラインメントを合わせるリファレンスのオブジェクト  
                     select = "backbone",  # 計算対象を主鎖にする  
                     ref_frame = 0    # リファレンスのフレームインデックス  
                        )
RMSD_analysis.run() 

run()で解析を実行すると、結果はresults.rmsdに「[[frame, time(ps), RMSD(A)], [...], [...]]」という情報を含むnumpy.ndarrayの形で保存されています。

print("type : ", type(RMSD_analysis.results.rmsd))
print("1st data : ", RMSD_analysis.results.rmsd[0])
print("total data number : ", len(RMSD_analysis.results.rmsd))

# type :  <class 'numpy.ndarray'>
# 1st data :  [0.0000000e+00 0.0000000e+00 6.0600393e-07]
# total data number :  50

matplotlibでプロットしましょう。0 psはリファレンスそのものなので除きます。

time_ps = RMSD_analysis.results.rmsd[1:, 1]
backbone_RMSDs = RMSD_analysis.results.rmsd[1:, 2]

plt.plot(time_ps, backbone_RMSDs)

plt.title("change of the RMSD")
plt.xlabel("Time (ps)")
plt.ylabel("RMSD [$\AA$]")
plt.show()

右肩上がりのプロットになりました。600 ps以降は1.1Å前後をうろうろしている感じです。

ついでにRMSDの分布もプロットしてみましょう。seabornのkdeplotを使ってカーネル密度推定法のプロットを行います。

import seaborn as sb

sb.kdeplot(RMSD_analysis.results.rmsd[1:, 2], shade=True)

plt.title("distribution of the RMSD")
plt.xlabel("RMSD [$\AA$]")
plt.show()

トラジェクトリ初期のリファレンスとのRMSDが小さい状態は短く、後半のRMSDが1.1Å前後まで広がった状態が多くを占める様子がわかりました。

4-4. RMSF

4-4-1. RMSFって?

つづいてRMSF(平均二乗揺らぎ, Root-mean-square fluctuation)を計算します。これは分子内原子の平均位置からの分散を示す量で、平衡状態からのトラジェクトリ構造の揺らぎを表すそうです。*7

4-4-2. RMSFの計算

RMSFの計算にはMDAnalysis.analysis.rms.RMSF()を使います。引数には計算したいatomgroupを指定します。

アラインメントが取られたタンパク質を使うことが想定された機能なので先のalign.AlignTrajを実行していない場合は処理を忘れずに!

Cα原子についてRMSFを計算する場合はこんな感じです。

# Cαの選択
C_alphas = md_universe.select_atoms("name CA")

# RMSF計算
RMSF_analysis = rms.RMSF(C_alphas)
RMSF_analysis.run()

run()で解析を実行すると、結果はresults.rmsfnumpy.ndarrayとして格納されます。RMSDのときと同じですね。

結果をプロットしましょう。Cα原子の残基番号はresnumsで取り出せます。

# 残基番号の取り出し(resnums)
residue_numbers = C_alphas.resnums

# プロット
plt.plot(residue_numbers, RMSF_analysis.results.rmsf)

plt.title("C_alpha RMSF")
plt.xlabel("residue number")
plt.show()

両末端でRMSF(揺らぎ)が大きくなっているのは納得として、850番 ~ 900番にも揺らぎのピークがあります。具体的にはEGFRタンパクのどのあたりの構造なのでしょうか?面白そうなのでもう少し調べます。

4-4-3. 文献計算例との比較

今回のシミュレーションのRMSF解析の結果が妥当か?文献計算例と比較してみましょう。

まず、具体的なタンパク質構造に戻って確認してみます。PyMolで該当の領域(850 ~ 900)を色分けしてみました。

だいたいこの辺。

キナーゼ構造でいうところのActivation loop(A-loop)に相当しそうです。

(ググって最初に出てきた文献と)結果を比較してみましょう。

Tamirat MZ, Koivu M, Elenius K, Johnson MS
Structural characterization of EGFR exon 19 deletion mutation using molecular dynamics simulation.
PLoS ONE 2019(14): e0222814

EGFRの構造を漫画で描くと以下で、DFGモチーフ(D855, F856, G857)に続いてA-loopと呼ばれる構造があるそうです。

上記論文は野生型(Wild EGFR)と変異型(欠損の変異、Δ746ELREA750)についてMDシミュレーションを行なって、変異が構造に与える影響を調べたもののようです。詳細は論文をご参照いただくとして、RMSFの評価を引用します。

野生型、変異型ともにA-loopでRMSFが大きくなっています!また、値としても1〜 2Åということで近い値でした。

素人目には類似の結果が得られているようにみえるのですが、専門家の目からみるとどうなんでしょうか?ご教示いただけると嬉しいです。

4-5. 回転半径(radius of gyration)

最後に回転半径の計算を行ってみます。こういう感じ指標です。*8

回転半径の算出にはAtomGroupに紐づいたメソッドであるradius_of_gyration()を使えば良さそうです。例えばCα原子の回転半径を求めるならこんな感じ。

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

# Cαの選択
C_alphas = md_universe.select_atoms("name CA")

# 回転半径の計算
C_alphas.radius_of_gyration()

# 19.713873743719788

19.7Åと出ました。select_atoms()で目的の原子群を選択してAtomGroupを作っておいてメソッドを適用するのはこれまで何度か出てきた通りですね。

トラジェクトリ全体を通した変化を見たい場合は、Universeオブジェクトのtrajectory属性についてforループを回してフレームイテレーションするのが簡単です。*9*10

# フレームを格納するリスト
frame_list= []
# 回転半径を格納するリスト
rgyr_list = []

# トラジェクトリのフレームごとに回転半径を計算してリストに格納
for ts in md_universe.trajectory:
    f = ts.frame
    rgyr = C_alphas.radius_of_gyration()
    
    frame_list.append(f)
    rgyr_list.append(rgyr)

# プロット
plt.plot(frame_list, rgyr_list)

plt.title("change of the radius of gyration")
plt.xlabel("frame")
plt.ylabel("rgyr [$\AA$]")
plt.show()

「1 ns ( = 0.002(ps) x 500,000 (steps))」のシミュレーションが計50フレームなので「1フレーム = 20 (ps) = 10,000 (steps)」です。上下はあるものの途中までは継時的に増加し、19.8Å付近で平衡にあるように見えます。RMSDで見た挙動に似ていますね。

回転半径の分布もプロットしてみましょう。

sb.kdeplot(rgyr_list, shade=True)

plt.title("distribution of the radius of gyration")
plt.xlabel("rgyr [$\AA$]")
plt.show()

回転半径の分布はRMSDの分布と比較して分布の肩の部分が高くなっています。同じトラジェクトリでも指標によってばらつきの度合いが変わるのは面白いですね。(・・・生物学的な解釈はわかりませんが、、、)

尚、ドキュメント(MDAnalysis User Guide:Radius of gyration)には回転半径を求める関数を自作する例が書かれています。ご興味のある方はご参照ください。

5. おわりに

以上、今回の記事ではMDシミュレーショントラジェクトリの解析を試みました。

冒頭に挙げたコロキウム資料を参考に、①Jupyter notebook上でNGLViewを使った動画の確認を行ったのち、②MDAnalysisを使ったRMSD、③RMSF、④回転半径の算出を行いました。

解析操作で忘れがちな注意点としてはアラインメントをとることとトラジェクトリのフレーム位置のリセットでしょうか?

何はともあれ、素人でもチュートリアルで動かせるOpenMMMDAnalysisはとても素晴らしいですね!オープンソースコミュニティに感謝です!

個人的にはRMSFでEGFRキナーゼA-loopの揺らぎが確認できた点が非常に面白かったです。また、RMSDと回転半径で分布が変わってくる点も興味深かったです。

今回取り上げたもの以外にもMDトラジェクトリの解析には色々とあるようです。続きはCM②のあとで!

よく理解せずにまま実行動かしているので間違いが多いと思います。ご指摘いただければ幸いです。

ではでは

*1:下図のMDAnalysisのロゴはドキュメントから。イラストはいらすとやさんより。

*2:NGLviewのDevelopment versionのインストール方法の部分です。

*3: フラッグ--sys-prefixは現在の(アクティベートされてる)仮想環境にkernelをインストールするためのもののようです。

*4:参考:How to properly use --sys-prefix。kernelのインストールがnglviewのアクティベートに必須かはよくわかりませんでした。

*5:引数「in_memory=True」はmobileで指定したトラジェクトリを永久にメモリ上のトラジェクトリ(in-memory trajectory)に変えるもののようです。そのためか、一度align.AlignTrajを実行したUniverseオブジェクトに再度実行してもファイルの書き出しが行われませんでした。そこでここでは再度オブジェクトを作り直しています。

*6:図の記載は理化学研究所 計算科学研究センターのeラーニング「計算生命科学の基礎6 生体系分子シミュレーションの新展開」における横浜市立大学 池口満徳先生(研究室HP)の講義動画(YouTube)とスライド資料(SlideShare)などを参考にしています

*7:こちらも図の記載は先のeラーニング資料などを参考にしています。

*8:こちらも図の記載は先のeラーニング資料などを参考にしています。

*9:参考:MDAnalysis:Overview over MDAnalysis

*10:フレームのイテレーションについてはMDAnalysis UserGuide:Slicing trajectoriesなども参考になります。

OpenMMをステップバイステップで 〜Part 5 シミュレーションプログラムの中身を理解しよう〜

OpenMMの使い方を順番にたどる記事のPart 5です。前回までで、GUIOpenMM Setupを利用して作成したファイル(①タンパク質構造、②シミュレーション実行プログラム)を用いて、Google Colab上でシミュレーションが実行できることを確認できました。

今回は「シミュレーションプログラムの中身がどのようになっているか?」、 OpenMMのコードの書き方について眺めて見たいと思います。

また、コードの解釈に先立って「OpenMM プログラムの全体像(各クラス間の相関関係)」を把握したいと思います。

1. プログラムの全体像(相関図)

1-1. Diagram of classes in OpenMM 6.0

OpenMMのプログラムの全体像について参考になる図がありました。バージョン 6.0*1についてのものですが、クラス構成は以下のようになっているそうです。

OpenMMのクラス構成

上図はSimTKプロジェクトにあるOpenMMDocumentsで公開されているOpenMM Class Diagramです。*2

各クラスの図形の意味は左下の点線枠内にあります。

OpenMMはPythonスクリプトで実行できるようになっていますが、内部の仕組みはC++で書かれているものもあります。これらは「API Layer Class」として緑色ひし形でダイアグラムに示されています。SystemContextStateといったクラスです(OpenMM C++ API)。

一方で、Pythonで書かれているものは「App Layer Class」 として黄色円形で示されています。SimulationReportersなどです(OpenMM Python API)。

また、「App Layer Class」のうちファイルを読み込んで働くものは灰色円形となっています。PDBFileや他のMDソフトで作成した入力(AmberPrmtopFIle etc.)などです。

さらに、各クラスに紐づいたメソッドや生成されるデータは、それぞれ赤枠カード黄枠カードになっています。

これらのクラスの関係性を矢印でつないでいるのが上図のダイアグラムです。*3

1-2. ちょっと簡略化

ダイアグラムがかなり盛りだくさんなので少し簡略化してみます。

メインになりそうな部分を白抜き表示にしています。それぞれこういう役割のクラスのようです。

ドキュメントに基づくクラスの説明

説明文だけだとよくわかりません・・・。

2. OpenMMのクラスとMD計算上の役割?

2-1. イメージ

文章だとクラスの役割がわかりにくいので、MDシミュレーションの実行工程とムリやり対応づけてみましょう。

まず、MD計算の流れはこんな感じでした。

これをふまえて、クラスと対応づけるとこんな感じでしょうか???

定義した計算System)でシミュレーションSimulation)を行って、得られた系の時間発展(状態の"文脈", Context)から知りたい時点の情報State)を適宜取り出すイメージです。

2-2. くどく説明

もう少しくどく書きます。

① まず、シミュレーションを行う系の定義をSystemクラスで行います。「どんなタンパク質?水分子はいくつ?どんな箱に置くの?力場は?…」みたいな基本的な条件設定のイメージです。

② ついで、計算の実行に関する設定項目(「時間発展の積分計算はどうするか?」)や出力方法(「結果のデータやトラジェクトリはどうやって出力する?」)に関する設定を、先のSystem)と結びつけるためにSimulationクラスを利用します。

③ このSimulationを使ってエネルギー最小化minimizeEnergy())やトラジェクトリの生成(プロダクション、step())を行うことになります。

④ また、MD計算は時間変化を追うシミュレーションなので、系の状態のデータを格納するためにContextクラスが必要になります。例えばシミュレーション開始時点ではタンパク質の各原子の位置情報(pdb.positions)を設定したりします。

⑤ その上で、特定の時点での状態の情報を取り出すためには別途Stateクラスを利用することになります。シミュレーションの要である情報を格納したContext変数には直接アクセスできないStateを経由してアクセスする)仕組みになっているようです。

⑥ 最後に、シミュレーション全体の時間変化やデータはreporterを使って出力することになるようです。

3. サンプルコードを眺めよう!

OpenMMの概観がわかったところで、OpenMM HPにあるサンプルコードを例にPythonの書き方を眺めてみましょう!

3-1. パラメータの復習

コードをみる前に、過去記事でとりあげたOpenMM SetupGUIツール)で設定したパラメータのイメージを貼っておきます。

設定項目のイメージ

たくさんあって大変です。。。

3-2. サンプルコード

サンプルコードではどうなっているでしょうか?

from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout

pdb = PDBFile('input.pdb')
forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter('output.pdb', 1000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))
simulation.step(10000)

(設定項目の多さから考えれば)10数行でMD計算の実行ができるというシンプルさがOpenMMの良さですね!

3-3. コードの中身

3-3-1. ステップ分割

ではコードの中身を見てみましょう。import部分は省略して、残りは以下のように分割できそうです。

サンプルコードを工程で分割

青い点線で「⓪系の準備、①エネルギー最小化、②シミュレーション本番」の3つのステップに区切っています。ここでは「平衡化」に相当する操作がないことにご留意ください。

各コードの右側には「今まで見てきたOpenMMのどのクラスに相当するか?」を書いています。

SystemSimulation、そしてContextと順番に構築して「⓪ 系の準備」を行います。ついで「② エネルギー最小化(minimizeEnergy())」を経て「③ シミュレーション本番(step(10000))」を行っています。

データの出力設定(reporters)は「保存したいステップ(③ シミュレーション本番)の直前」に行うというのもポイントになりそうです。

なお、サンプルコードでは状態の確認などを行っていないのでStateクラスが出てこないことにご留意ください。

3-3-2. 系の準備-1

各ステップを細かく見ていきましょう。まずは系の準備です。

最初に、Systemオブジェクトをタンパク質やその他粒子のトポロジーに基づいてForceFieldから構築します。

まず、(水分子の配置などの前処理済みの)PDBファイルの読み込みPDBFile())と、力場のパラメータを定義したXMLファイルの読み込みForceField())を行います。

ここで、ForceFieldは水モデルの力場も合わせて読み込んでいることにご注意ください(amber99sb.xmltip3p.xml)*4

読み込んだらトポロジーpdb.topology)と、(相互作用)の計算方法を引数にとってcreateSystem()することでSystemオブジェクトが構築できます。

なお、力の計算の設定は非結合相互作用(粒子・メッシュ・エバルト法、nonbondedMethod=PME)やカットオフ距離拘束条件などです。

3-3-3. 系の準備-2

次に、Simulationオブジェクトを構築します。これはtopologysystemintegratorを引数にとることで構築できます。

イラストはいらすとやさんより

Systemの構築に使ったtopologyがもう一度出てくるのが少しわかりづらいです。

積分計算方法(Integrator)としてサンプルではLangevinIntegratorを利用しています。これはランジュバン動力学を利用するもので、温度制御ができます。

LangevinIntegratorの引数は(temperature, frictionCoeff, stepSize)で、サンプルでは熱浴の温度を表す引数に「300*kelvin」が設定されています。

なお、nanometerkelvinといった単位が利用できるのは最初に以下を実行しているからです。

from openmm.unit import *

単位系が使えるのは便利なので忘れずに実行したいですね。

3-3-4. 系の準備-3

系の準備の最後はsimulationからのcontextの構築です。

位置についてContextを作成(イラストはいらすとやさんより)

ここではPDBファイルの情報のうち「位置 pdb.position」を利用します。

以上で系の準備は終わりです。

3-3-5. エネルギー最小化

ついで、系のエネルギーを最小化します。

minimzeEnergy()simluationに行うだけです。シンプル!

3-3-6. シミュレーション本番

最後にシミュレーション本番です。トラジェクトリのプロダクションです。

Reporterを追加して計算のstepを踏んでいく(イラストはいらすとやさんより)

step()により、指定したタイムステップのシミュレーションが実行できます。

ただし、ここではトラジェクトリを出力したいので、時間発展させる前に種々のreporterオブジェクトをreporters.append()によりsimulationに追加しています。

出力内容としてPDBReporter()は指定したファイルパスに、一連のフレームを(一定間隔ごとに)PDBファイルとして書き出します。

また、StateDataReporter()はエネルギーや温度といった状態の情報を書き出します。

それぞれの細かな設定はマニュアルをご参照ください。

以上がサンプルコードの中身でした。

4. OpenMM Setupで作成したコードを眺めよう!

4-1. サンプルコードとの違いの概略

サンプルコードを見ることで、シミュレーション実行のためのコードの流れがつかめました。

もう少し長いコード例として、過去記事OpenMM Setupで作成したものを眺めてみましょう!

コード全体の流れは、予めオブジェクトの作成やメソッドのための色々な引数変数に代入する形で定義しておき、後ほどシミュレーション実行コードを書くという順番です。

また、サンプルコードのシミュレーションとの大きな違いとして以下の点が挙げられます。

  • ① NPTアンサンブルで圧力制御も行う
  • ② シミュレーション本番前に平衡化のステップも行う

細かな違いとしては以下の点が挙げられます。

  • ③ Reporterの種類の違い
  • ④ 計算機環境(platform)の設定を明示的に行っている

4-2. コード全体を眺めよう

違いをふまえた上でコード全体を眺めましょう!

コードの途中にコメントアウトする形で説明を追記しました。

# 必要なライブラリのインポート

from openmm import *
from openmm.app import *
from openmm.unit import *

# ファイルの読み込み(PDBと力場XML → Systemのベースになるファイル)

pdb = PDBFile('3poz-processed.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

# Systemの構築に使うパラメータ変数の設定  
## 非結合相互作用の計算方法や拘束条件の定義

nonbondedMethod = PME
nonbondedCutoff = 1.0*nanometers
ewaldErrorTolerance = 0.0005
constraints = HBonds
rigidWater = True
constraintTolerance = 0.000001

# 積分計算(Integrator)の設定に必要なパラメータ変数の設定  
## ここではLangevinMiddleIntegrator
dt = 0.002*picoseconds
temperature = 300*kelvin
friction = 1.0/picosecond

# アンサンブルNPTとしているため圧力制御も行う(サンプルコードとの違い)
## 圧力制御のパラメータ変数を設定
pressure = 1.0*atmospheres
barostatInterval = 25

# シミュレーションを行うための設定
## シミュレーション時間の設定  
## 平衡化も行うため別途タイムステップを定義している(サンプルコードとの違い)
steps = 500000
equilibrationSteps = 50000

## 計算環境をPlatformで設定している(サンプルコードとの違い)
platform = Platform.getPlatformByName('CUDA')
platformProperties = {'Precision': 'single'}

## 予めreporterとして何を使うかも変数に定義している
dcdReporter = DCDReporter('trajectory.dcd', 10000)
dataReporter = StateDataReporter('log.txt', 1000, totalSteps=steps,
    step=True, speed=True, progress=True, potentialEnergy=True, temperature=True, separator='\t')
checkpointReporter = CheckpointReporter('checkpoint.chk', 10000)

# シミュレーションの準備  
## ForceFiledからSystemを構築する
print('Building system...')
topology = pdb.topology
positions = pdb.positions
system = forcefield.createSystem(topology, nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
    constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance)

## 圧力制御の機能をaddForceする(サンプルコードとの違い)
system.addForce(MonteCarloBarostat(pressure, temperature, barostatInterval))

## 積分計算の設定
integrator = LangevinMiddleIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)

## simulationの構築
simulation = Simulation(topology, system, integrator, platform, platformProperties)
simulation.context.setPositions(positions)

# XML結果出力ための設定
with open("system.xml", mode="w") as file:
    file.write(XmlSerializer.serialize(system))
with open("integrator.xml", mode="w") as file:
    file.write(XmlSerializer.serialize(integrator))

# エネルギー最小化の実行
print('Performing energy minimization...')
simulation.minimizeEnergy()

# 平衡化のステップ(サンプルコードとの違い)
print('Equilibrating...')
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(equilibrationSteps)

# シミュレーション本番の実行の前にreporterオブジェクトを追加する

print('Simulating...')
simulation.reporters.append(dcdReporter)
simulation.reporters.append(dataReporter)
simulation.reporters.append(checkpointReporter)
simulation.currentStep = 0

# シミュレーション本番の実行
simulation.step(steps)

# シミュレーション最後の状態をXMLで書き出す
simulation.saveState("final_state.xml")

ずいぶん長くなっているように見えますが、実際の操作としてはサンプルコードとそれほど変わっていなさそうです。

4-3. コードの違いを少し詳しく

折角なのでサンプルコードと違う部分を少し詳しくみてみましょう。

4-3-1. NPTアンサンブル(圧力制御)

まず、アンサンブルの違いです。

NPTアンサンブル達成のためMonteCarloBarostatという圧力制御を系に加えています。これは、Systemを構築後にaddForce()の形で追加することで利用可能です。

圧力制御の方法(イラストはいらすとやさんより)

温度制御に関わるIntegratorSimulationを構築する際の引数であったのと異なり、その前の「系の構築段階で力として加える」という点は要注意です。

4-3-2. プラットフォーム

続いて計算機環境の設定です。Platformクラスを利用して設定しています。

計算機プラットフォームの設定(イラストはいらすとやさんより)

ここではGoogle ColabのGPUCUDA)を利用するためgetPlatformByName('CUDA')となっています。

Platformはそのプロパティ{'Precision': 'single'})と合わせてSimulationを構築する時に引数として与えられています。

Simulationが計算の実行や出力のさまざまなオブジェクトを取りまとめる役割だったことを考えると納得できそうです。

4-3-3. 平衡化

最後に平衡化ステップを加えることによる違いです。

平衡化の計算自体はシミュレーション本番と同様に、simluationに対して必要な時間に相当するstep()を行えばOKです。

初速分布を与えて平衡化(イラストはいらすとやさんより)

ただし上図のように、OpenMM Setupで作成したコードではさらにsetVelocitiesToTemerature()というものをcontextに与えています。

これは「系の全粒子に対して、設定の温度ボルツマン分布に従うようにランダムな速度を与える」という役割を担います。エネルギー最小化後の系に初速を与えたと考えればよさそうです*5

contextが速度を含む系の状態を保持するクラスだったことをふまえれば、このような追加方法になるのも納得できそうです。

また、平衡化のstep()の計算が終わったのち、シミュレーション本番のstep()に入る前に以下のコードを実行しています。

simulation.currentStep = 0

平衡化の状態は保持しつつ「ステップのインデックスcurrentStep)をゼロ」にすることで、シミュレーション本番のトラジェクトリ出力(reporter)に合わせていると考えられそうです。

以上がコードにおける主な違いでした。

5. クラスの振り返り

サンプルコード、およびOpenMM Setupで作成したコードの振り返りが終わりました。

記事全体のまとめも兼ねてコードにでてきたOpenMMクラスのダイアグラムを作ってみましょう。

こんな感じでしょうか?

(すこしアップデートした)簡略化クラスダイアグラム

一番最初に引用したOpenMM Class Diagramの半分くらいは抑えられたかな?、という印象です。

設定項目とクラスの種類が多いので関係性がわかるとコードを書きやすくなりそうですね。

6. (おまけ)OpenMMのアーキテクチャ

最後に脱線ですが、OpenMMの構成(アーキテクチャ)の図を引用しておきます*6

OpenMMのアーキテクチャ

OpenMMは計算を実行するためのコードの部分では、計算機の環境(プラットフォーム)を意識しなくて良いように設計されているそうです。

利用者はアーキテクチャ上部のPublic Interface(Public API)の使い方を把握しておけばよく、これはプラットフォームと独立なコードになっています。このPublic APIが上で見たSystemContextなどです*7

プラットフォームに依存した部分は、アーキテクチャー下部のPlatform Abstraction Layer(Low Level API)やComputational Kernelsです。最下層が実際に計算を実行するハードウェアに関する部分だそうです。

独自の計算手法を開発する場合などは下層の理解が必要となるそうですが、そもそもそんな上級者はこんなブログ読みませんよね。。。

以上のように、プラットフォーム依存の部分独立の部分分かれた構成となっているので、(貧弱な)自分のパソコン上でもGoogle ColabのパワフルなGPUを使った計算のためのプログラム作成が簡単にできたのでしょうか???

よくわかりませんがアーキテクチャって言ってみたかっただけです。ごめんなさい。

7. おわりに

以上、今回はOpenMMでシミュレーションを実行するためのコードの書き方について、「① (OpenMM HPの)サンプルコード」と「② OpenMM Setupで作成したコード」を例に眺めてみました。

また、コードの中身を眺めるに先立って、OpenMM Class Diagramを利用して各クラス間の相関関係を大まかに把握し、OpenMMプログラム全体の概観の理解を試みました。

OpenMMは(ばりばりの計算の専門家でなくても)シミュレーションを実行できるように、Pythonで簡単にシミュレーションプログラムをかけるよう設計されているそうです。

実際に、MD計算に必要な設定項目の多さ多段階工程のややこしさからすれば、シンプルなコードを並べるだけで計算を実行できるようになっていました。

・・・とはいうものの、やはりある程度は複雑になってしまいます。

今回の記事で見たように、クラスダイアグラムのような形式で各操作の関係性がわかっていると動作するコードをちゃんとかけるようになるかもしれませんね?

以上、今回も間違いが多そうです。ご指摘いただければ幸いです。

ではでは!

*1: Class Diagram資料の公開日付は2014年のようです。2022年現在はOpenMM 7.7が最新版となっています。

*2:SimTKは医学生物学計算コミュニティのための無料のプロジェクトホスティングプラットホームだそうです(Whit is SimTK?)。

*3:Lee-Ping Wang博士によるワークショップ(YouTube)「Creating and Customizing Force Fields in OpenMM」での説明などを参考にしています。

*4:利用可能な力場のXMLファイルについてはOpenMM User Guide : 3.6.2. Force Fieldsなどを参照してください。

*5:サンプルコードで初速を設定していない理由はよくわからないです…。NPTアンサンブルのMonteCarloBarostatにも温度がある関係?

*6:OpenMM User Guide: 7.4. Architectural Overview

*7:OpenMM User Guide: 7.5. The OpenMM Public API

OpenMMをステップバイステップで 〜 Part 4-2 よくある質問(FAQ@GitHub)〜

OpenMMのシミュレーションを実行した前回の記事(part 4)で、結果の描画でタンパク質が境界から飛び出ているのが気持ち悪いと書きました。

OpennMMのGitHubページWikiFrequently Asked Questions」に全く同じ内容が書かれていました。全く気づいていませんでした。。。

よくある質問なら共有しておいた方が良いよね!ということで、DeepL翻訳(にして修正)したものを置いておきます。*1*2

f:id:magattaca:20220409170433p:plain
みんな同じところでつまずいてる!

1. 周期境界条件はどのように機能するの?

 非常によくある質問は、もとを探ると周期的境界条件の理解に関するものです。「なぜ水分子が周期境界条件の箱の外側に拡散するのか?」「なぜタンパク質が水のボックスの端からはみ出るのか?」といった質問がよくあります。 また、「特定の力が周期境界条件においてどのように適用されるか?(あるいは適用されないか?)」といったことに関する質問もあります。 一見、簡単な質問に見えますが、見た目よりもずっと複雑な問題なので、まず初めに説明しましょう。

 概念的に言えば、周期境界条件でシミュレーションを行う場合、すべての「粒子」は、空間全体で繰り返しパターンをとる無限の粒子のグリッドを表します。 周期的なボックスが10nmの立方体で、ある粒子が(1, 2, 3)にあるとすると、(11, 2, 3)、(21, 2, 3)、(-9, 12, 103)・・・にも粒子があることになります。これらの位置のうち、どの位置をリストアップするかは重要ではありません。 どれをとっても、全く同じ粒子の無限のグリッドを表しています。

 次に、理解すべき非常に重要なポイントがあります。

周期境界条件は、位置ではなく、変位量(displacement)に適用されます。

 もう一度10nmの立方体を例にとりましょう。(1, 0, 0) と (9, 0, 0) の2つの粒子間の変位量はどうなるでしょうか? 周期境界条件がない場合、変位量は(8, 0, 0)となります、周期境界条件がある場合、(-2, 0, 0)となります。 もし、2番目の粒子の位置が(-1, 0, 0)や(99, 0, 0)であったとしても、それは問題ではありません。変位量は同じになります。 周期境界条件は、粒子間の変位量に依存する力を計算するときに適用されます。

 もう一つ、非常に重要なポイントがあります。

周期境界条件は、非結合相互作用のみに適用され、結合相互作用には適用されません。

 2つの粒子間の結合を考えてみましょう。それぞれの粒子は実際には無限の粒子のグリッドを表していますが、だからといって、最初の粒子のすべてのコピーが2番目の粒子のすべてのコピーに結合しているわけではありません。 最初の粒子の各コピーは、2番目の粒子の1つの特定のコピーと結合しており、それ以外のものは結合していません。 もし第2の粒子の異なるコピーが、結合しているコピーよりも何らかの方法で近づくことができたとしても、結合が突然「切れて」異なるコピーと結合をつくりなおすはずだということにはならないのです。

 (実際には、本当に必要であれば、結合相互作用に周期境界条件を適用することは可能です。結合相互作用に対して setUsesPeriodicBoundaryConditions() を呼び出すだけでできます。 ですが、通常は行わないことで、特別な状況でのみ適切です。 例えば、一方の端がもう一方の端の次の周期的なコピーに結合しているような無限分子をシミュレートしたいといった場合があげられます。)

 これは拘束条件や、実際には結合相互作用の一種となっている非結合相互作用の例外にも適用されます。非結合相互作用は通常、互いに結合している粒子間(1-2、1-3、1-4相互作用)では除外されるか、縮小されます。 これは直接結合している粒子にのみ適用され、その粒子の他の周期的なコピーには適用されません。

 その結果、2つの異なる周期的なコピーの間で分子を壊してはいけないということになります。先に、粒子が(9, 0, 0)にあるか(-1, 0, 0)にあるかは重要でないと書きました。しかし、もしその粒子が(1, 0, 0)にある粒子と結合しているならば、それは重要なことなのです。 一方は結合長が8であることを示し、もう一方は結合長が2であることを示します。ほとんどのMDパッケージはこれと同じ要件を満たしているので、通常は問題ありません。 しかし、時々、すべての原子が独立して1つの周期ボックスにまとめられているファイルに出くわすことがあり、分子が真っ二つになっていることがあります。 MDTrajのmake_molecules_whole()関数など、このようなファイルを修正するためのツールがあります。

 ここまではシミュレーションの実行に関することでしたが、次に出力ファイルの保存について説明します。 MDエンジンにとって、(0, 0, 0)を中心とする分子は(100, -30, -80)を中心とする分子と全く同一かもしれませんが、トラジェクトリを可視化しようとする場合、異なる分子が多くの異なる周期的な箱に散らばっていると何が起こっているか分からなくなってしまいます。 このため、Context.getState() には、返されたすべての粒子の位置を一つの周期的なボックスにまとめるオプションがあります。 同様に、レポーター関数(reporter)は周期系のトラジェクトリを書き出すときに自動的にこれを行います。ボックスにまとめる作業は常に分子全体に対して行われます。すべての分子の中心は1つのボックスの中にありますが、その分子のいくつかの原子はその外側にはみ出していることがあります。これが、タンパク質の一部だけが水に囲まれているように見えることがよくある理由です。

 次のような重要なポイントがあるため、可視化はさらに複雑になっています。

一意に定義された「周期的なボックス 」は存在しない。

 例えば、ボックスの幅が10である場合、粒子の位置はすべて0と10の間にあるべきということでしょうか? それとも-5から5の間でしょうか?どちらの方法も様々なプログラムで共通する慣習になっています。 タンパク質が原点を中心とする初期構造を読み込んで、すべての分子を0から10の間に置くようなトラジェクトリを書きだすと、タンパク質のほとんどが水の外に出てしまい、反対側の水のボックスに大きな穴が開いているように見えます。 これは非常によくある問題なので、多くの解析プログラムには特にこの問題に対処する機能が備わっています。例えば、MDTraj の image_molecules() 関数や Cpptraj の autoimage コマンドなどがあります。これらのツールは、すべての粒子を一つの箱に入れて、溶質を中心に、その周りに水を配置するような賢い方法をとります。

2.「粒子の位置がNaNです」「エネルギーがNaNです」といったエラーが出るのはなぜ?

 NaNとは「not a number」の略で、「数字ではない」という意味です。NaNは、ゼロで割ったり、負の数の平方根を取ったりするような、実数値がうまく定義できない計算の結果をコンピュータが表現する方法です。このエラーが発生した場合、一般的にはシミュレーションが破綻(blown up)していることを意味します。このようなエラーは様々な理由で発生しますが、一般的には力や速度が大きすぎることが原因です。その結果、粒子が急激にジャンプし、次のタイムステップでさらに大きな力が発生し、数ステップですべての粒子の位置がNaNに変わってしまいます。

 シミュレーションが破綻する原因には、さまざまな種類の問題があります。より一般的なものをいくつか紹介します。

  • タイムステップが大きすぎる場合。 タイムステップが大きすぎる場合、シミュレーションはしばらくの間、正常に実行されたように見えますが、その後、破綻します。タイムステップが大きすぎる場合、すぐに不安定になります。

  • 初期コンフォメーションがエネルギー最小化されていない場合。 例えば、結晶構造を初期コンフォメーションとして使用している場合、座標が力場の最小値になっていない可能性があり、最初のタイムステップで非常に大きな力を発生させます。これは通常、局所的なエネルギー最小化を実行することで修正することができます。

  • 周期的なボックスのベクトルが正しく設定されていない場合。 例えば、Amberのprmtopとinpcrdファイルから開始する場合によくある間違いです。どちらのファイルも周期的ボックスのベクトルを指定しており、お互いに異なる場合があります。粒子座標があるボックスサイズを想定していても、シミュレーションが異なるボックスサイズを使用している場合、粒子間の激しい衝突を引き起こす可能性があります。

  • 拘束条件が正しく定義されていない場合。 拘束条件の適用には行列を逆行列にする作業が含まれます(明示的な場合もあれば、暗黙的な場合もあります)。この行列の条件が悪いと、シミュレーションが破綻する原因になります。単純な原因としては、同じ粒子ペアに2つの異なる拘束条件が課されている場合が挙げられます。また、同時にすべてを満たすことができない矛盾した拘束条件がある場合にも起こりえます。その他、より複雑な拘束条件のセットも、条件付けが悪い行列につながる可能性があります。

  • 力の定義が正しくない場合。 これは、系(System)を直接構築しているとき、あるいは系に余分な力を追加しているときに、最も起こりやすい問題です。適用する力が大きすぎる場合、シミュレーションが破綻する可能性があります。また、カスタマイズした力を作成する際には、直接NaNを生成するような式を書かないように注意してください(負の値を取りうるの値の平方根を取る式など)。

 このようなエラーは様々な原因で発生するため、デバッグが困難な場合があります。 最初に確認することは、シミュレーションがすぐに(最初の数ステップで)破綻するか、はじめは長時間実行されるかです。問題が初期コンフォメーション(最小化されていないなど)または系の定義(拘束条件や力など)のいずれかにある場合、通常すぐに破綻します。局所的なエネルギー最小化を実行してみて、それで解決するかどうかみてください。 一方で、ステップサイズが大きすぎる場合、はじめはしばらく実行されることがあります。その場合は、まずステップサイズを小さくしてみてください。

 エラーが断続的なものではなく再現性がある場合、通常は何が起こっているのかを正確に追跡することができます。 まず、非常に大きな力がかかっている粒子を探すことから始めてください。そして、その力がどこから来ているのかを特定します。その力は、結合性相互作用によるものか、それとも非結合性相互作用によるものでしょうか? 具体的な相互作用を特定できますか? 例えば、2つの粒子が近すぎるのでは?

 これらのエラーをデバッグする別のアプローチは、エラーの原因となるものを系統的に削除することです。力を1つずつ取り除いてみてください。 すべての拘束条件を取り除いてみる。あるいは、周期境界条件を無効にしてみる。少しの辛抱と忍耐で、大抵の場合、原因を特定することができます。

3. 「残基のテンプレートが見つかりません」というエラーが出るのはなぜ?

このエラーは、ForceField上で createSystem() を呼び出したときに発生することがあります。このエラーの意味を理解するためには、OpenMMで力場がどのように動作するかを少し知っておく必要があります。

 力場を使って系(System)を構築する最初のステップは、すべての原子のタイプを特定することです。 原子のタイプは非常に固有のものである場合があります。例えば、全ての炭素原子が同じであるわけではありません。力場は、α炭素と芳香族側鎖の炭素とでは、おそらく異なるパラメータを使用するでしょう。 原子タイプは、残基をテンプレートにマッチングさせることで決定されます。力場は、認識できるすべてのタイプの残基のためのテンプレート定義を含んでいます。createSystem()を呼び出すと、渡されたトポロジー(Topology)のすべての残基をループして、それぞれをテンプレートにマッチングさせようとし、そのテンプレートで指定された原子タイプを割り当てます。

 残基は2つの情報、すなわち各原子の元素とそれらをつなぐ結合(残基の原子と別の残基の原子をつなぐ結合も含む)の情報に基づいてテンプレートにマッチングされます。残基の名前など他の情報については気にしません。(ただし、後述するように、誤った残基名は問題を引き起こす可能性があります。) 元素と結合だけが重要です。

 残基にマッチングするテンプレートを見つけることを妨げる多くの問題があります。ここでは、より一般的なものをいくつか紹介します。

  • 力場のXMLファイルが省略されている。 これは最も単純な可能性の1つです。ForceFieldオブジェクトを作成し、どのファイルを読み込むかを指示したときに、単に1つだけ抜けている可能性があります。例えば、Amber14とCHARMM36の定義では、水を他の力場とは別のファイルに入れ、異なる水モデルを簡単に選択できるようにしています。そのファイルを指定しなかった場合、力場はトポロジー(Topology)の水分子に一致するテンプレートを持ちません。

  • 非標準の残基がある。 AmberやCHARMMのような標準的な力場は、タンパク質や核酸のような様々な標準的な分子に対応しています。しかし、 Topologyにリガンドやその他の低分子が含まれている場合、 力場にはそのためのパラメータがない(つまり、 テンプレートがない)可能性があります。

  • トポロジーに原子がない。 例えば、 結晶構造には通常、水素原子が含まれておらず、 またフレキシブルな領域では重原子が欠落していることがよくあります。 これらの原子がない場合、トポロジーは力場のテンプレートとマッチングしません。 OpenMM-Setup は、この種の問題の多くを自動的に修正することができます。

  • 水モデルが Topology と矛盾している。 4サイトおよび5サイトの水モデルは、各水分子が仮想サイトのための「余分な粒子」を含んでいることを想定しされています。これがない場合(あるいは逆に、あるのに3サイトの水モデルを使おうとした場合)、テンプレートがマッチングしません。Modeller.addExtraParticles()を使って、この問題を修正することができます。

  • 鎖の終端が正しくない。 タンパク質、そして核酸はさらに、末端の構造が多種多様になっていることがあります。鎖の末端の残基は、中間位置の残基とは違っていて、異なるテンプレートが必要です。多くの場合、力場は各鎖の終端の構造について1つか2つしかサポートしていません。もしTopologyの末端が異なるものである場合、テンプレートにマッチングしません。

  • 非標準のPDBファイルを持っている。 PDBは標準化されたフォーマットであり、specificationで正確に定義されています。 残念ながら、多くのPDBファイルはその仕様に従っておらず、様々な問題を引き起こしています。

 最後のポイントについては議論が必要です。PDBファイルには様々な非標準の可能性があるからです。テンプレートは原子の元素と原子間の結合という2つの要素に基づいてマッチングされることを思い出してください。これらの一つでも不正確に再構築されると、テンプレートのマッチングでエラーが生じます。

 PDBファイルの各ATOMレコードには、元素を指定するフィールドがあります。残念ながら、そのフィールドが空白のままのファイルもあります。その場合、PDBFileクラスは原子の名前から元素を推測しようとします。通常は正しく推測することができますが、常にそうとは限りません。

 結合はより複雑です。PDB の仕様では、標準的な残基と異質原子(heterogen)を区別しています。標準的な残基(アミノ酸ヌクレオチド、水)には結合が指定されていません。ファイルを読み込むソフトウェアが自動的に結合を追加することになっています。その他の分子はCONECTレコードで結合が明示的に指定されているはずです。これらのレコードがない場合、PDBを読むソフトウェアはどの原子が結合しているのかわからないので、ForceFieldはテンプレートをそれらにマッチングさせることができません。

 標準残基には独自の問題があります。それらは名前によって認識されます。 PDB Chemical Components Dictionary ではすべての残基と原子の標準的な名前が指定されています。これにより、例えば、ALAという残基はCAという原子とCBという原子の間に結合があるはずだということがPDBを読むソフトウェアにはわかるわけです。残念ながら、ファイルによっては残基や原子の名前に非標準的なものを使っている場合があります。そのような場合、PDBFileはどこに結合を追加すればよいのかわかりません。

 これが、原子や残基の名前が重要になる理由です。ForceFieldはテンプレートをマッチングするときに名前を使いません。しかし、ファイルが非標準的な名前を持っていると、PDBを読むソフトウェアは正しい結合を追加することができず、ForceFieldがテンプレートにマッチングすることができなくなります。

 PDBFileは非標準的な名前に対応するように努めています。よく使われる非標準的な名前を認識するための大きなテーブルが用意されています例えば、水分子の標準的な名前は HOH ですが、それ以外にも WATSOLTP3 などの名前も受け付けます。同様に、水分子の酸素の標準的な名前はOですが、OWOH2という名前も受け付けます。これらの名前は、過去に様々なプログラムが書いたPDB ファイルで誤って使用したもので、OpenMM はそれらを探すことができます。しかし、OpenMMが認識できない他の非標準的な名前を使用している場合、結合は正しく追加されず、ForceField はテンプレートをマッチングさせることができなくなります。

4. SystemやForceに加えた変更が無視されるのはなぜ?

 一旦Contextを作成すると、そのContextが作成されたSystemはイミュータブル(不変のオブジェクト)として扱われます。そのシステムに対して更に変更を加えても無視されます。これには、ForcesやVirtualSitesなど、Systemに含まれる他のオブジェクトも含まれます。

 Contextの作成にはコストがかかります。GPU カーネルコンパイル、スプラインの表形式関数へのフィッティング、拘束条件行列の構築と逆行列の計算など、シミュレーションの開始時に行う必要がある先行作業が多くあります。その過程で、Systemからの情報は多くの異なる場所にコピーされ、多くの場合、高度に処理された形でコピーされます。SystemまたはForceに変更が加えられると、多くの内部データ構造が無効になる可能性があります。そのため、OpenMM では、効率と混乱を避けるために、既存のContextはSystemに対するすべての変更を無視するというルールになっています。Contextを作成すると、必要なものはすべてSystemからコピーされ、それ以降は参照されません。

 このルールは、状況によっては妨げになる場合があります。シミュレーションの過程で何かを変更することが本当に必要な場合があります。このような状況に対応するために、スピードと柔軟性のトレードオフの関係にある複数の方法が提供されています。

 最初の方法はグローバルパラメータです。グローバルパラメータの値はSystemではなくContextに保存されるので、setParameter()を呼び出すことでいつでも変更することができます。カスタマイズしたフォースは全て、グローバルパラメータを定義することができます。 パラメータに依存する形でエネルギーを定義した場合、パラメータの新しい値を設定するだけで、いつでも簡単にエネルギーを変更することができます。 NonbondedForceでは、パラメータのオフセットで同様のことを行うことができます。グローバルパラメータを定義することで、粒子のセットを一度に協調的に変更することができます。

 次に、ほとんどのフォースではupdateParametersInContext()メソッドが用意されています。これはSystemの限られた部分を更新するための制御された方法を提供します。例えば、粒子のセットの電荷を変更することができます。これを行うには、NonbondedForceに対してsetParticleParameters()を呼び出します。 他の変更と同様に、これは既存のContextでは無視されます。 しかし、次にupdateParametersInContext() を呼び出すと、変更されたパラメータ値が既存のContextにコピーされます。これは、グローバルなパラメータを変更するだけの方法よりもはるかに遅いですが、完全に新しいContextを作成するよりもはるかに速いです。この方法では、Forceのすべての条件を変更できるわけではないことに注意しましょう。updateParametersInContext() メソッドが何を更新し、何を更新しないかについては、各メソッドの API ドキュメントを参照してください。

 最後の選択肢は、Context上でreinitialize()を呼び出すことです。 これにより、内部のデータ構造をすべて破棄し、一から作り直します。 これは、事実上まったく新しい Context を作成するのと同じことです。非常に遅いですが、絶対に全てを変更することができます。

 内部データ構造がすべて破棄されるため、位置、速度、パラメータ値などの状態の情報も失われます。引数 reinitialize(preserveState=True)を渡すことで、それらの情報を保持するように要求することができます。これを実行すると、まずメモリ上にチェックポイントを作成し、次に Context を再初期化(reinitialize)し、最後にチェックポイントを再読み込みします。ほとんどの変更に対して、これは簡単に状態の情報を保持する方法を提供しますが、もしチェックポイントを有効でなくするような方法でSystemが変更された場合は、失敗して例外を発生させます。例えば、System内の粒子の数が変化した場合、チェックポイントをロードすることができなくなり、例外が発生し、状態の情報が失われます。

以上が、2022年4月現在のFAQです。

*1:OpenMMのGitHubにはライセンスファイルが見当たりませんでしたが、調べた限りOpenMMがMITあるいはGPLのようなので翻訳させていただいています

*2:イラストはいらすとやさんより

OpenMMをステップバイステップで 〜 Part 4 ColabでOpenMM シミュレーションを実行 〜

OpenMMの使い方を順番にたどる記事のPart 4です。過去の記事で、タンパク質構造の前処理(Part 1Part 2)とシミュレーション実行のためのスクリプトPart 3)を手に入れるところまでが終わりました。

あとは計算を実行するだけ!でも大変そうだからColabでやってしまおう!ってなお話です。*1

f:id:magattaca:20220403195515p:plain
オープンソースとColabで低コストMD

1. 準備

1-1. 使用ファイル

今回使用するファイルは以下の2つです。

  • PDB ID: 3POZをMDシミュレーション用に処理したファイル(3poz-processed.pdb)
  • ② OpenMMシミュレーションを実行するスクリプトrun_openmm_simulation.py

それぞれの作成手順は過去の記事をご覧いただければと思います。

処理済みのPDBファイルはこんな見た目です。

f:id:magattaca:20220403195536p:plain
計算する系の見た目

1-2. スクリプトの修正

Google Colabでシミュレーションを実行する前に、前回作成したスクリプトを少し修正しました。

修正内容は2点です。

  • ① シミュレーションの途中経過(残り時間 etc.)をColabノートブック上に表示させる
  • ② シミュレーション後の最終状態をPDBx/mmCIFフォーマットでも出力する

1-2-1. 途中経過の出力

MDシミュレーションのトラジェクトリ生成段階(プロダクション)は計算時間がかかるので、「計算そのものが進んでいるか?」「どれくらい時間がかかるのか?」心配になります。

計算のログ自体はログファイル(log.txt)に書き出すようにスクリプトで設定されていますが、計算途中でも確認できた方が安心できそうです。

そこでColabノートブックのセルの出力にも同様の内容を書き出すことにしました。

f:id:magattaca:20220403195734p:plain
計算の残り時間がわかった方が嬉しい

上図のように、プロダクションを行う「simulation.step(steps)」の直前にもう一つStateDataReporter()を追加しました。標準出力(stdout)に進捗情報を表示するためのコードです。

from sys import stdout
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, progress=True, remainingTime=True, speed=True, totalSteps=steps, separator="\t"))

1-2-2. 最終状態のPDBx/mmCIF出力

前回のスクリプト作成では、シミュレーション後の出力ファイルをXMLにしていました。

# Write file with final simulation state

simulation.saveState("final_state.xml")

最終状態の系を可視化するならPyMolで読み込めた方が楽だよね!と、いうことで構造を出力するようにしました。

f:id:magattaca:20220403195837p:plain
PDBx/mmCIFフォーマットで結果を出力

OpenMM SetupOutputタブで「Save final simulation statePDBx/mmCIF (no velocities)」とすると上記のようになります。

この状態のスクリプトを使っても良いですし、シミュレーション終了後に同じコードを実行してもファイルを書き出すことができます。

state = simulation.context.getState(getPositions=True, enforcePeriodicBox=system.usesPeriodicBoundaryConditions())
with open("final_state.pdbx", mode="w") as file:
    PDBxFile.writeFile(simulation.topology, state.getPositions(), file)

2. Google Colabでシミュレーション

2-1. 手順

それではGoogle Colabに計算してもらいましょう!

手順は以下です。*2

f:id:magattaca:20220403195920p:plain

くどくど書いてますが、インストールしてスクリプト(.py)を実行するだけです。

2-2. 実行

では実際に計算を進めていきましょう!

2-2-1. GPUに設定

Colabにログインして新しいノートブックをつくったら忘れずにハードウェアをGPUにしておきます。

f:id:magattaca:20220403200150p:plain
いつも途中で思い出すやつ

以下をセルで実行して、GPUに切り替わったことを確認することもできます。

!nvidia-smi

f:id:magattaca:20220403200220p:plain

2-2-2. OpenMMのインストール

次に、OpenMMをインストールします。

OpenMMはcondaがあれば簡単にインストールできるので、最初にminicondaをインストールします。

!wget -qnc https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
!bash Miniconda3-latest-Linux-x86_64.sh -bfp /usr/local
!rm Miniconda3-latest-Linux-x86_64.sh

wgetLinux用のminicondaを取ってきて、「/usr/local」にインストールしています。終わったら不要なファイルをrmで削除してます。

次に、OpenMMをconda-forgeチャンネルからインストールします。

!conda install -c conda-forge openmm=7.6 python=3.7 -y 

インストールできたらOpenMMが使えるようにパスを通しておきます。

import sys
sys.path.append('/usr/local/lib/python3.7/site-packages/')

忘れててModuleNotFoundErrorで怒られました。

f:id:magattaca:20220403200320p:plain
パスを通さないと怒られる

2-2-3. ファイルアップロード

Colabに計算で使用するPDBファイル3poz-processed.pdb)とスクリプトrun_openmm_simulation.py)をアップロードします。

f:id:magattaca:20220403200354p:plain

Colabの左側でクリックするだけ。前からこんなに簡単でしたっけ?

2-2-4. スクリプトの実行

あとはスクリプトを実行するだけ!

%run run_openmm_simulation.py

標準出力に経過を出すように書き換えたので進捗が確認できます。

f:id:magattaca:20220403200425p:plain
計算の進行度がわかる

今回はトラジェクトリの長さを1nsとしていますが、プロダクションだけでおおよそ40分程度かかりました(平衡化などを含めると小一時間くらい)。

2-2-4. 結果のダウンロード

計算が終わると左側ファイルに色々と新しいファイルができています。

f:id:magattaca:20220403200450p:plain
出力いろいろ

google.colabfilesを使って、簡単にローカルPCにダウンロードできます。*3

from google.colab import files
files.download("trajectory.dcd")
files.download("log.txt")
files.download("checkpoint.chk")
files.download("system.xml")
files.download("final_state.xml")
files.download("final_state.pdbx")

Colabはランタイムが切断されるとファイルが全て消えてしまうので、とりあえず目につくものをダウンロードしておきました。

容量が大きいファイルはブラウザが「本当にダウンロードして良いの?」と聞いてくる場合があるので許可してください。*4

PCではなくGoogle Driveに保存したい時はマウントしてコピーすれば良いそうです。

例えば以下でGooge Driveのマイドライブディレクトリに「trajectory.dcd」が保存されます。

from google.colab import drive
drive.mount('/content/drive')

!cp trajectory.dcd /content/drive/MyDrive/trajectory.dcd

Colabの左側に出ているデフォルトのディレクトリ名はcontentとなっているようです。上記でcontent/driveGoogle Driveがマウントされます。

パスに困った時は左側に表示されるディレクトリで「パスをコピー」すればOKです。

f:id:magattaca:20220403200531p:plain
階層構造がよくわからない

3. PyMolで可視化

結果のファイルが手に入りました!取り急ぎPyMolで可視化してみましょう。

シミュレーションの最終時点のファイルは「final_state.pdbx」です。OpenMM Setupのデフォルトのファイル名では拡張子が「pdbx」となっているので「.cif」に名前を変えてからPyMolに入れました。

f:id:magattaca:20220403200613p:plain
final stateが変

境界からタンパク質がはみ出したぞ????なんでや???

よくわからないですが動画にしました。構造を読んだ状態でトラジェクトリ(trajectory.dcd)をPyMolに読み込ませれば良さそうです。

f:id:magattaca:20220403200718g:plain

かわいい!!・・・けど、はみ出してる。。。

4. Center of Mass Motion?

境界からタンパク質が出てしまっている原因がよくわからないので検索するとCenter of Mass Motionという用語に行き当たりました。OpenMM User Guide3.6.12 Removing Center of Mass Motionという項目があります。

曰く、デフォルトではOpenMMアプリケーションは「CMMotionRemover」を系に付与しており、「各ステップ毎に重心の動きすべて取り除くことで、時間発展に伴い系全体がドリフトしてしまうことを防ぐ」のだそう、、、。

ドリフトが何かよく分かりませんが、タンパク質が境界外に出るのはかなり「drift漂流)」感のある動きです。

上記項目では、「removeCMMotion=Falseを設定すればドリフトを許容することもできるけど、かなりレアなシミュレーションだよねー」と書いてありました。

デフォルトから変えていないので重心は動かない設定のはずなのですが・・・。心当たりのある方、ご教示いただけると嬉しいです。

5. おわりに

以上、今回はGoogle Colab上でOpenMMの計算を実行する方法について見てみました。手順がくどくて読みにくい記事になってしまいましたが、記憶力の悪い私の備忘録ということでお許しください。

得られた結果を可視化するとタンパク質が境界外に出ているというなんとも気持ちが悪いものでした。何はともあれ、一応は動きのある結果が得られたということで、Colab上での計算手順は確認できたのでヨシ!

・・・困った。

以上、今回もかなり間違いが多そうです。お詳しい方ご教示いただけると嬉しいです。

ではでは!

6. 追記

境界からタンパク質がはみ出ることについてよくわからないと書きましたが、OpenMMGitHub Frequently Asked Questionsにそのものズバリ書いてありました。

Why does my protein stick out past the edge of the water box?

「問題ないよ!周期境界条件(periodic boundary conditions)をきちんと理解していないだけだよ!」ってな感じみたいです。

シミュレーション本番においても、計算後のトラジェクトリ描画においても、「境界をどのようにとるか?」「分子をどのように描画するか?」によって見え方が異なってくるようです。

GitHubのIssue #2688 Simple example of using SystemBuilder to simulate protein-ligand complexの議論中でも、タンパク質が境界から飛び出た図が引用されていますが、特におかしい挙動ではない(The trajectory produced looks just fine)と書かれています。

以上、同様の描画が得られて気になる方はご参照いただければと思います。

*1:イラストはいらすとやさんより

*2:Making it rainの手順を参考にさせていただいています。

*3:参考:Google Colabのファイルのアップロードとダウンロード

*4:ダウンロードが始まらないと思ってたら、ポップアップが出てました。気づかなかった。。。

OpenMMをステップバイステップで 〜 Part 3:GUIでシミュレーション条件の設定 〜

OpenMMの使い方を順番にたどる記事のPart 3です。Part 1Part 2の記事でMD計算のためのタンパク質構造の前処理が完了しました。今回の記事ではシミュレーション条件本体のセットアップの流れを追ってみたいと思います。

利用するツールはOpenMMGUIOpenMM Setupです。なお、今回の記事では設定のみで計算そのものは実行しません。*1

1. OpenMM Setupでできること

まず最初にOpenMM Setupでできることの概要です。

前々回の記事でPDBファイルを読み込んで構造の前処理の方法を確認しました。この時はPDBの出力で一度終わりとしましたが、そのまま続けてシミュレーション条件を指定するスクリプトrun_openmm_simulation.py)の作成保存、さらに実行Run Simulation)までを行うことができます。

画面構成はこんな感じ。

f:id:magattaca:20220402183941p:plain
シミュレーション条件設定の画面構成

上図のように、左側で選んだ条件が右側のスクリプト逐次的に反映されて設定ファイルが出来上がっていきます。

前処理と同じく指示に従って条件を選択していくだけ!これなら私でもできそう!

2. MD計算の流れと設定項目のおさらい

シミュレーション条件を設定していく前に、MD計算の流れと設定項目のおさらいをしておきましょう。

計算の流れはざっと以下です。

f:id:magattaca:20220402184012p:plain
MD計算のおおまかな流れ

準備直後の(⓪)は水分子やイオンを適当に配置しただけなので、最初にエネルギー最小化(①)します。ついで各原子にランダムに速度を与え、徐々に系全体を設定温度まで引き上げていきます(平衡化(②))。最後に平衡状態の系においてアンサンブルを達成するようにトラジェクトリを生成します(シミュレーション本番(③))。*2

設定の必要な項目は大体こんな感じ。*3

f:id:magattaca:20220402184033p:plain
各工程で設定が必要となる項目

ごちゃごちゃしてますが、系の初期状態と平衡状態(のアンサンブル)を定めて、そこに至る計算方法分子力場インテグレータ)や、状態の制御方法サーモスタットバロスタット)を設定します。あとは希望のシミュレーション時間(ステップ数)計算すれば良さそうです。

また、MD計算はヘビーなので計算環境GPUをつかえるか?)やログ、結果(トラジェクトリ)の出力管理といった項目も大事な設定となってくるようです。

これからOpenMM Setup4つのタブで順番に設定を確認していきますが、おおよそ以下の対応になっています。

f:id:magattaca:20220402184052p:plain
各タブと設定内容の対応

では順番に見ていきましょう!

3. OpenMM Setupでシミュレーション条件を設定しよう!

3-1. 注意!~ 入力手順とデフォルト値の違い

設定を始めていく前に少し注意点があります。OpenMM Setupでは入力の手順で少し表示内容が変わります。

具体的にはPDBファイルを入力とした時、「① 未処理の構造を入力としてOpenMM SetupでClean upしたあと、続けてシミュレーションの設定に進んだ場合(Part 1 記事でご紹介した手順)」と、「② すでに処理済みPDBファイルを入力として直接シミュレーションの設定に進んだ場合」です。

f:id:magattaca:20220402184128p:plain
手順によって初期パラメータが異なる

上図左のように前者(①:Yes. Let's clean it up now.)では、シミュレーションの設定項目のほとんどが空欄でひとつずつ自分で選んでいきます。一方、後者(②:No. my file is all ready to simulate.)ではデフォルトの設定(数値)が最初から記入されていて、変更したい箇所だけを変えればOKです。

Clean upをした後でデフォルトの数値設定が知りたい場合は、処理後のPDBファイルを一度保存して、再度読み込んで最初からはじめればOKです。

3-2. Systemタブ

それではひとつずつタブを見ていきましょう。まずはSystemタブです。

「System(系)」とついていますが、ここでは主に長距離相互作用静電相互作用)の計算方法を設定します。

先の項(3-1.)で見た通り、タンパク質構造の前処理からそのまま進んだ時は項目の大半が空欄です。Systemでは以下のようになっています。

f:id:magattaca:20220402184201p:plain
Systemタブ(初期値がほぼ空白の場合)

Nonbonded MethodConstraintsといった項目があります。右側のスクリプトでは「# System Configuration 以下」に対応していることがわかります。

左側の各欄にカーソルを合わせると「何を設定する項目か?」ヘルプが表示されます。

f:id:magattaca:20220402184223p:plain

今回はTeachOpenCADDT019 · Molecular dynamics simulationやデフォルト値を参考に以下のように設定してみました。

f:id:magattaca:20220402184245p:plain
Systemの条件設定例

まず、① 長距離静電相互作用(nonbondedMethod)はPME(粒子・メッシュ・エバルト法、エバルトの方法 - Wikipedia)です。*4

カットオフ距離nonbondedCutoff)はPME法のうち、「粒子部分(実空間における短距離ポテンシャルの直接和)」に関するパラメータで、「直接計算する距離をどの程度の長さにするか?」を定めるカットオフ値のことだと思います。*5

Ewald法の誤差許容範囲ewaldErrorTolerance)は、「Ewald法で和を取る際に、切り捨て計算することで生じる力の端数の誤差」に、おおよそ等しくなるように設定します。デフォルトでは0.0005です。*6

拘束条件constraints)では、結合長や結合角に関する拘束をオプションで設定することができます。HBondsとした場合は、「水素原子を含む全ての結合の結合長」を拘束します。他にはAllBondsHAnglesといった選択肢があります。*7

拘束条件の誤差許容範囲constraintTolerance)は、「拘束条件について誤差をどの程度許容するか?」で、1x10-6としています。*8

なお、スクリプト側にある「⑥ 剛体としての水ridigdWater)」は水分子の扱いを決めるパラメータです。OpenMMでは水分子を結合長・結合角ともに拘束された、完全に剛直なもの(rigid)として扱います。Falseにするとflexible waterとなります。

以上がSystemタブでした。*9

3-3. Integratorタブ

つづいてIntegratorタブです。積分計算やアンサンブルについて設定します。

こんな感じで設定してみました。

f:id:magattaca:20220402184313p:plain
Integratorタブの設定例

まず、① 積分計算の時間間隔Step Size, dt)は「0.002ピコ秒(2フェムト秒)」としました。

統計的アンサンブルStatistical Ensemble)には、「Constant pressure, temperature, NPT」と「Constant volume, temperature, NVT」の選択肢があります。今回はNPTとしました。

尚、GUIのOpenMM Setupでは「アンサンブルを達成するための計算方法(温度や圧力を一定に保つための計算方法)」を選択する項目は無いようです。NPTを選択するとデフォルトで「温度LangevinMiddleIntegrator圧力MonteCarloBarostat」が設定され、スクリプトに反映されます。

温度temperature)は「300K(27℃)」にしました。

摩擦係数Friction Coefficient, friction)は温度制御に関わるパラメータです。デフォルトのLangevinMiddleIntegratorはランジュバン動力学に基づくものなので、ランジュバン方程式の摩擦項の係数をパラメータとして設定します。*10*11

圧力pressure)の単位は「atm(標準大気圧)」です。

圧力調整の間隔Barostat Interval, barostatInterval)は、「何ステップごとに圧力調整を行うか?」です。①で「Step size : 0.002 ps」としているので、⑥で間隔を「25 steps」とすると「0.002 x 25 = 0.050 ps」毎に調整を行います。*12*13

以上がIntegratorタブでした。

なお、上図のスクリプト部分は変数の設定のみです。計算方法の設定は以下のようになっていました。

f:id:magattaca:20220402184339p:plain
積分法、圧調整法のスクリプト

積分LangevinMiddleIntegratorや圧調整が MonteCarloBarostat となっていることが確認できますね。

3-4. Simulationタブ

続いてSimulationタブです。シミュレーション時間や計算環境の設定を行います。

Google ColabのGPUを利用するつもりで、以下のように設定してみました。

f:id:magattaca:20220402184357p:plain
Simulationタブの設定例

シミュレーションの長さはトラジェクトリを生成する時間(production)で、単位はstepです。

Integratorタブで積分Step Sizeを「0.002 ps」としました。従って、上図のように「500,000 steps」とすると、「0.002 x 500,000 = 1,000 ps = 1 ns」となり、1 ナノ秒のシミュレーションとなります。

平衡化の長さも同様にstepで長さを指定します。上図では「50,000 steps」(0.1 ナノ秒)としました。

尚、上記の時間設定はかなり適当なので参考にしないでください。目的にあった適切な値を選択してください。*14

プラットフォームでは計算環境を設定します。選択肢にはReferenceCPUCUDAOpenCLがあります。Google ColabではGPUを使用できるのでCUDAとしています。*15

精度はプラットフォームでCUDAあるいはOpenCLを選択した場合に選べるパラメータで、「singlemixeddouble」があります。

精度よりも速度を重視するなら「single(単精度)」で、逆に精度を重視するなら「double(倍精度)」です。「mixed」はをsingle、積分をdoubleで計算する混合モードです。*16

以上がSimulationタブでした。

3-5. Outputタブ

最後にOutputタブです。結果の出力を設定します。色々な形式で出力できて、項目数が多いので分割して見ていきます。

f:id:magattaca:20220402184549p:plain
Outpuタブの設定例(1)

OpenMMではトラジェクトリをPDBPDBx/mmCIFあるいはDCDフォーマットで保存することができます。*17GUIのOpenMM Setupでは、このうちDCDフォーマットのみを選択できます。*18

まず、① ボックスにチェックを入れて「トラジェクトリを保存する」を選択し、ファイル名と「どれくらいの間隔(step)ごとに結果を出力するか?」を指定します。

上図では「trajectory.dcd」ファイルに「10000 steps (20 ピコ秒)」毎に状態を書き込むことになります。

また、計算のログファイルを残すこともできます(②)。プロダクションタイムのトラジェクトリ以外のさまざまなデータを保存することができます。

保存する内容はData to Writeの中から選べ、エネルギーや速度、温度などがあります。上記では、「log.txt」ファイルに「1000 steps(2 ピコ秒)」毎に書き込みます。*19

①、②で設定した内容はスクリプトではそれぞれ①DCDReporter( )、②StateDataReporter( )に引数として反映されています。

f:id:magattaca:20220402184612p:plain
Outpuタブの設定例(2)

OpenMMではシミュレーションの進行具合をチェックポイントファイルというバイナリファイル(.chk)で残すことができます(③)。サーバーエラーなどでシミュレーションが途中で止まってしまった時に、チェックポイントファイルを読み込んで計算を再開することができるそうです。

上図では、「checkpoint.chk」ファイルに「10000 steps (20 ピコ秒)」毎に上書き保存しています。スクリプトでは③CheckpointReporter( )です。*20

f:id:magattaca:20220402184632p:plain
Outpuタブの設定例(3)

また、シミュレーションの状態(位置、速度 etc.)をXMLファイルに書き出すこともできます。このXMLファイルを再度読み込むことでシミュレーションを続けることができます。*21

chkには「同じ計算環境、OpenMMのバージョンでしかシミュレーションを再開できない」という欠点がありますが、XMLにはその欠点がありません。一方で、XMLにはサイズが大きくなるという欠点があります。

上図では、シミュレーションの設定(④)と、最後の状態(⑤)をそれぞれXMLで保存しています。

スクリプトでは、④シミュレーションの設定「system.xmlintegrator.xml」はテキストファイルに書き込む形式で書かれており、⑤最後の「final_state.xml」はsaveState( )というメソッドを使った形式になっています。

なお、⑤はXML以外に、chkPDBx/mmCIFフォーマットを選択することができます。

以上がOutputタブでした。

3-6. スクリプトを保存

これで4つのタブでの設定がすべて終わりました。最後に作成されたスクリプトを保存しておきましょう。

f:id:magattaca:20220402184658p:plain
スクリプトGet!

Save Scriptをクリックすると「run_openmm_simulation.py」というファイルが保存されます。

4. おわりに

以上、今回はOpenMM Setupを使ったシミュレーション条件の設定をみてみました。GUIでクリックしていくことで逐次的にスクリプトが生成され、同時に確認できるというのはとても分かりやすくて良いですね!

以前の記事と合わせて、タンパク質構造とシミュレーションスクリプトが揃ったことになります。あとはこれが本当に動くのかどうかが心配です・・・。 早速試したいところですが、記事が長くなってしまったので続きは次回に。

記事の途中、設定の説明を試みていますがOpenMM User Guideを適当に訳しただけなので、意味のわからないところ、間違っているところが多いと思います。ご指摘いただけると嬉しいです。

おしまい。

*1:今回の記事はOpenMM User Guide3. Runnning Simulationの記述を参考にしています。

*2:参考:MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜

*3:イラストはいらすとやさんより

*4:プルダウンで選べるPME以外の選択肢は図の右上

*5:OpenMM User Guide 3.1. A First Exampleでは「cutoff for direct space interactions」となっています。湯どうふさんのYouTubeの解説動画「分子シミュレーション実行時のクーロンポテンシャルの効率的な計算法:Ewald法とParticle Mesh Ewald法【理論化学、計算科学】」がとても分かりやすいです。

*6:OpenMM User Guide 3.6.5. Nonbonded Interactions

*7:OpenMM User Guide 3.6.6. Constraints

*8:OpneMMのGitHub Constraint Tolerance Guidelinesでは「1x10-5と1x10-6のどちらが良いか?」ディスカッションされています。「温度一定なら1x10-5、エネルギー一定なら1x10-6」といった意見が出ています。

*9:今回は設定しませんでしたがHydogen Mass Repartitioningは「重原子に結合する水素原子の重さを変えて、動きの速い水素原子を遅くする」ものだそうです。OpenMM User Guide 3.6.7. Heacy Hydrogens

*10:ランジュバン動力学 - wikipedia

*11:OpenMM User Guide 3.6.8. Integrators

*12:OpenMM User Guide 3.6.10. Pressure Coupling

*13:MDの圧力調整におけるモンテカルロ法(Monte Carlo barostat)については、森野キートスさんのブログ記事「AmberでMD計算するときの圧力制御法の違い」などが参考になります。

*14:OpenMM Setupのデフォルト値では、「Step Size :0.004 ps、Simulation Length :1000000 steps、Equilibration Length:1000 steps」のようです。

*15:それぞれの選択肢の詳細はOpenMM User Guide 7.7. Platformsをご参照ください。Referenceはパフォーマンスは考えずに読みやすいコードにしたもので、他の選択肢はそれぞれのプラットフォームに合わせてパフォーマンスを最適化しています。

*16:各プラットフォームで選択できる精度を含む個別のパラメータについてはOpenMM User Guide 10. Platform-Specific Propertiesをご参照ください。

*17:OpenMM User Guide 3.6.13. Writing Trajectories

*18:DCDはMD計算のトラジェクトリでよく使われるフォーマットのようです。参考:DCD Plugin, Version 1.11

*19:OpenMM User Guide 3.6.14. Recording Other Data

*20:OpenMM User Guide 3.6.15. Saving Simulation Progress and Results

*21:OpenMM User Guide 3.6.15. Saving Simulation Progress and Results

molzipで色々なリンカーのPROTACをつくろう!

RDKitのとても便利な機能molzipというのを知りました。複数のフラグメントをまとめてひとつ分子を構築できるというもので、① @dr_greg_landrum先生の記事(R-group decomposition and molzip)と、② @iwatobipen先生の記事(Build peptide from monomer library from ChEMBL)に使い方が詳しく紹介されています。

①の記事ではRグループ分析と組み合わせた分子ライブラリの構築方法を、②の記事ではモノマーを組み合わせたペプチドライブラリの構築方法を紹介してくださっています。

ユニットを組み合わせて作るタイプの分子に威力を発揮しそうですね!ユニットといえば他に思いつくのは・・・PROTAC!!

というわけで、molzipを使って色々なリンカーのPROTACを作ってみたいと思います!また、構造データの取得元としてPROTAC-DBを参照します。

f:id:magattaca:20220327214104p:plain

※ 注)この記事でやることは基本的に上記2記事と一緒です。遊んでみたかっただけです。手抜きでごめんなさい。

1. molzipの基本

1-1. 基本的な使い方-1

まずはmolzipの基本的な使い方を確認しましょう!molzipのドキュメントはこちら(RDKit documentationmolzip)。

こんな感じで別々の分子(Molオブジェクト)を1つにまとめる関数です。

f:id:magattaca:20220327203926p:plain

「どこで分子どうしをつなぐか?」はオプションの引数(MolZipParams)で設定されます。デフォルトではMolzipLabel.AtomMapNumberとなっています。

あらかじめそれぞれの分子の結合を作りたい原子目印(ダミーアトム + ナンバリング)をつけておいて、それらの原子同士をつなぐ(ジッパーでとめる)感じです。

・・・チャックするってもう言わないの???

1-2. molzipの使用例-1

簡単な使用例を見てみましょう!Buchwald-Hartwigクロスカップリングをしたいと思います。

f:id:magattaca:20220327204012p:plain

C-N結合をつくりたいので、上図のように反応点をダミーアトム(*)にしてナンバリングしたMolオブジェクトを作りましょう!

from rdkit import Chem

mol1 = Chem.MolFromSmiles('c1ccccc1[*:1]')
mol2 = Chem.MolFromSmiles('C1CCCCN1[*:1]')

Chem.Draw.MolsToImage([mol1, mol2])

f:id:magattaca:20220327204054p:plain

できました!ワイルドカードアスタリスク(*)にコロン(:)を挟んで数字を置いて、角格好([ ])で囲えばOKです。*1

あとはmolzipするだけ。

product = Chem.molzip(mol1, mol2)
Chem.Draw.MolToImage(product)

f:id:magattaca:20220327204129p:plain

マッピングした原子同士で結ばれた分子ができました!

1-3. 基本的な使い方-2

さて、2つのMolオブジェクトをつなぐという使い方でした。

molzipには他にも「複数のフラグメント1つのMolオブジェクトに入ったものを対象にして分子を作る」、といった使い方もあります。

例えば、以下のようにSMILESをピリオド(.)を挟んで繋げれば複数のフラグメントが入ったオブジェクトを作れます。

mol3 = Chem.MolFromSmiles('[*:1]c1ccccc1.C1N([*:1])CC(C)(C)N([*:2])C1.c1nc([*:2])ccc1')
Chem.Draw.MolToImage(mol3)

f:id:magattaca:20220327204154p:plain

ベンゼンピリジンが近接していて見づらいですが、フラグメント3つ2種類のアトムマッピング(1,2)を含んでいます。

これをmolzipするとこんな感じ。

product2 = Chem.molzip(mol3)
Chem.Draw.MolToImage(product2)

f:id:magattaca:20220327204222p:plain

それぞれ指定したマッピングの箇所でつながった分子ができました!3つ以上のフラグメントをつないで分子を構築したい場合に便利そうです。

2. この記事でやること

molzipの使い方がわかったところで、この記事の目標です。

タンパク質分解誘導キメラ分子(Proteolysis targeting chimeda、PROTAC)は、① 分解したい対象のタンパク質に結合するWarheadと、② ユビキチンリガーゼに結合するE3 ligand、そしてをつなぐLinker(③)の3つの要素から構成されています。

f:id:magattaca:20220327204315p:plain

何を分解するか?(Warhead)」「どうやって分解するか?(E3 ligand)」という点に注目がいきがちですが、細胞内で分解機能を発揮するためには適切なリンカーを見つけ出すことが重要だそうです。

上図に引用した下記文献では、それぞれの標的ごとに試行錯誤しなければいけないリンカーについて焦点を当てて議論されています。(オープンアクセス, CC BY 4.0*2

Troup RI, Fallan C, Baud MGJ.
Current strategies for the design of PROTAC linkers: a critical review.
Explor Target Antitumor Ther. 2020;1:273-312.

例えばこんな構造があるそう。

f:id:magattaca:20220327204346p:plain

さらに組み合わせたり、長さが違ったり、色々あって大変だ!・・・それなら、どんな組み合わせ構造があるか全部作っちゃえばよくない?

というわけで、いろんなリンカーと残りのユニット(Warhead、E3 ligand)の組み合わせをmolzipでつくってみましょう!

3. PROTAC DBから構造を取得しよう

3-1. PROTAC DBって?

さて、いろいろな組み合わせ構造をつくろう!といっても、元になる構造群がないとどうにもなりません。

「どこかにまとまっていないかな?」・・・ということで素敵なデータベースがありました。その名もPROTAC-DB!!

文献は以下(オープンアクセス, CC BY 4.0

Gaoqi Weng, Chao Shen, Dongsheng Cao, Junbo Gao, Xiaowu Dong, Qiaojun He, Bo Yang, Dan Li, Jian Wu, Tingjun Hou,
PROTAC-DB: an online database of PROTACs,
Nucleic Acids Research, Volume 49, Issue D1, 8 January 2021, Pages D1381–D1387

Fig. 1を引用します。

f:id:magattaca:20220327204505p:plain

色々な文献からPROTACの例をあつめてきて活性データといっしょにまとめたデータベースのようです(PROTACにフォーカスしたChEMBLみたいな感じ?)。

さらに素晴らしいことに、それぞれのPROTACをWarheadLinkerE3 Ligandユニットに分割した構造のアノテーションもしてくださっているようです。これは便利!

3-2. PROTAC DBにアクセス!

早速アクセスしてみましょう。URLはこちら(「PROTAC-DB : http://cadd.zju.edu.cn/protacdb/ 」)

こんな感じでした。

f:id:magattaca:20220327204615p:plain

最新版のリリースは2021-04-17で、約2000個の分子が登録されています。

例としてArvinas社のBET degrader、ARV-771を検索してみましょう*3。結果のサマリーはこんな感じ。

f:id:magattaca:20220327204737p:plain

標的タンパク(BRD2, BRD3, BRD4)やE3ユビキチンリガーゼ(VHL)、種々の表記IUPAC、InChI etc.)といった情報があります。また、構造を3つのユニットに分けて分類していることも確認できました。

スクロールダウンすると、RDKitで計算したプロパティ(MW, logP, TPSA etc.)や、活性値といった情報も見ることができます。

f:id:magattaca:20220327204832p:plain

ウェスタンブロットの図まであるのにびっくり!元の文献にもReferenceをクリックで飛べます。素敵!

3-3. データをダウンロード

それでは目的の構造情報を取得しましょう!

下図のようにDownloadsをクリックして、移行先のページでDATA SOURCESから目的のユニットを選びます(今回はLinkers)。

f:id:magattaca:20220327204906p:plain

表示されるテーブルは標的タンパクごと(Compoud Group)にまとめられていて、SDF形式あるいはCSV形式でダウンロードすることができます。

今回はDB上にある全てのリンカーの情報が欲しいので「All行のlinker.csv」をダウンロードしました。

ダウンロードしたCSVファイルを確認してみます。

f:id:magattaca:20220327204947p:plain

上図のように構造(SMILES)、その他プロパティ計算値が格納されています。

SDFではなくCSVを選んだポイントは「Smiles_R列」です。結合箇所が「[R1][R2]」で示されていて分かりやすくなっています。SDFではこの情報がないため、「どこを結合箇所にすればよいか?」わからなくなっています。

先に見た通りmolzipの処理では、結合箇所の指定(AtomMapNumber)が大事でした。ですのでCSVファイルの「Smiles_R」のように結合箇所のわかりやすい情報の方が使いやすそうです。

何はともあれ、これでリンカーの情報が手に入りました。

4. molzipでいろいろなリンカーのPROTACをつくろう

4-1. 基本骨格の選択(ARV-771)

それでは色々なリンカーのPROTACをつくってみましょう!

まず、今回はLinker以外の①Warhead、②E3 ligandは固定したいと思います。

先にPROTAC-DBで見たARV-771を基本骨格として利用することにしましょう。

Raina K, Lu J, Qian Y, Altieri M, Gordon D, Rossi AM, Wang J, Chen X, Dong H, Siu K, Winkler JD, Crew AP, Crews CM, Coleman KG.
PROTAC-induced BET protein degradation as a therapy for castration-resistant prostate cancer.
Proc Natl Acad Sci U S A. 2016 Jun 28;113(26):7124-9.

PROTAC-DBの該当ページ(ID:328)からcanonical SMILESをコピペしてRDKitに読み込ませます。

ARV771 = Chem.MolFromSmiles('CC1=C(C2=CC=C([C@H](C)NC(=O)[C@@H]3C[C@@H](O)CN3C(=O)[C@@H](NC(=O)COCCCOCCNC(=O)C[C@@H]3N=C(C4=CC=C(Cl)C=C4)C4=C(SC(C)=C4C)N4C(C)=NN=C34)C(C)(C)C)C=C2)SC=N1')
Chem.Draw.MolToImage(ARV771)

f:id:magattaca:20220327205108p:plain

4-2. R-group分析で3つのユニットに分解

これを3つのユニットに分けましょう。リンカーをコア構造にしてR-group decompostionすれば良さそうです。参考のRDKit blog記事の通りやればOK!

まずコア構造をPROTAC-DBのARV-771 Linker構造をもとに作ります。

ARV771_linker = Chem.MolFromSmiles('NCCOCCCOC') 
Chem.Draw.MolToImage(ARV771_linker)

f:id:magattaca:20220327205154p:plain

R-group decompositionします(参考:RDKit ドキュメント rdRGroupDecompose)。

from rdkit.Chem import rdRGroupDecomposition

rgroups, _ = rdRGroupDecomposition.RGroupDecompose([ARV771_linker], [ARV771], asRows=False)

引数は①コア構造、②対象の分子の順番で、それぞれ複数の構造を入力できるようになっているので、リストにして渡しています。

結果は辞書型なので、keyを確認しましょう。

print(rgroups.keys())

# 出力
#  dict_keys(['Core', 'R1', 'R2'])

CoreR1R2とラベル付けされた3つに分解されました。

辞書のvalueはMolオブジェクトを格納したリストです。描画します。

Chem.Draw.MolsToGridImage([rgroups['Core'][0], rgroups['R1'][0], rgroups['R2'][0]], legends=['Core', 'R1', 'R2'])

f:id:magattaca:20220327205236p:plain

うまくできていそうです。

4-3. 結合部位の再ラベル化

このままリンカーとの連結に進みたいところですが、WarheadがR2、E3 ligandがR1となっています。PROTAC-DBと結合点のナンバリングがずれてしまうので、「Warhead : 1」「E3 ligand : 2」となるように修正しましょう。

WarheadのMolオブジェクトを確認すると、結合点のAtomSymbolR2になっています。

ARV771_warhead = rgroups['R2'][0]

for atom in ARV771_warhead.GetAtoms():
    if atom.GetSymbol() == "R2":
        for prop in atom.GetPropNames():
            print(prop, ':', atom.GetProp(prop))

# 出力
# molAtomMapNumber : 2
# dummyLabel : R2

上のように、プロパティとして「dummyLabel : R2」「molAtomMapNumber : 2」が付与されていることがわかります。SetPropを使ってラベルを「* : any atom」、ナンバリングを「1」にしましょう。

for atom in ARV771_warhead.GetAtoms():
    if atom.GetSymbol() == "R2":
        atom.SetProp('dummyLabel', '*')
        atom.SetProp('molAtomMapNumber', '1')

Chem.Draw.MolToImage(ARV771_warhead)

f:id:magattaca:20220327205348p:plain

できました!

同様にして、E3 ligandについてR1を「dummyLabel : *」「molAtomMapNumber : 2」に変更します。

ARV771_E3 = rgroups['R1'][0]

for atom in ARV771_E3.GetAtoms():
    if atom.GetSymbol() == "R1":
        atom.SetProp('dummyLabel', '*')
        atom.SetProp('molAtomMapNumber', '2')

Chem.Draw.MolToImage(ARV771_E3)

f:id:magattaca:20220327205440p:plain

見づらいですが、右端が「* : 2」になっています。

4-4. リンカーライブラリの作成

WarheadとE3 ligandが準備できたので色々なリンカー部位のデータベースを作成します。

PROTAC-DBからダウンロードした「linker.csv」をPandasのDataFrameに読み込みます。

import pandas as pd

linker_df = pd.read_csv('linker.csv')[['Compound ID', 'Smiles_R']]

f:id:magattaca:20220327205512p:plain

先に見た通り「Smiles_R列」に結合箇所を「[R1][R2]」で明記したSMILESが格納されています。これを「[* : 1]、[* : 2]」に置き換えましょう。

テキストを置換する関数をつくってmap()でDataFrameに適用します。

def linker_converter(smilesR):
    r1 = smilesR.replace('[R1]', '[*:1]')
    r2 = r1.replace('[R2]', '[*:2]')
    return r2

linker_df['smiles_dummy'] = linker_df['Smiles_R'].map(linker_converter)

f:id:magattaca:20220327205545p:plain

ダミーアトムに置換したSMILESが作成できていそうです。これでリンカーライブラリができました。

4-5. ユニットを連結しよう!

だいぶんと長くなってしまいましたが、やっと3つのユニットの基礎となる情報が揃いました。あとはひたすら繋いでいくだけです。

molzipの基本的な使い方で見た通り、3つのユニットをつなぐにはフラグメントをまとめて1つのMolオブジェクトにしておけばOKでした。

参考記事(R-group decomposition and molzip)では、RWMolオブジェクトに変更した後、Insert()してフラグメントをまとめていますが、ここでは簡単のためSMILESを連結する方式にします。

WarheadとE3 ligandをそれぞれSMILESに変換して、linkerのSMILESと「ピリオド(.)」を挟んでつないでいきます。

warhead_smiles = Chem.MolToSmiles(ARV771_warhead)
E3_smiles = Chem.MolToSmiles(ARV771_E3)

def linker_connector(linker_smiles):
    connected_smiles = warhead_smiles + '.' + linker_smiles + '.' + E3_smiles
    return connected_smiles

linker_df['protac_smiles'] = linker_df['smiles_dummy'].map(linker_connector)

print(len(linker_df))

# 出力
# 1099

f:id:magattaca:20220327210057p:plain

SMILESが長すぎて出来ているのかよくわからない。。。とりあえず1099行あるようです。

あとは1つずつ取り出してMolオブジェクトに変換した後、molzipしていきます。

molzipに失敗する場合もあるので、成功したものを「リスト:PROTACs」、失敗した場合のリンカーを「リスト:failed_linker」に格納していきます。

PROTACs = []
failed_linker = []
idx = 0

for smi in linker_df['protac_smiles']:
    fragments = Chem.MolFromSmiles(smi)
    try:
        protac = Chem.molzip(fragments)
        PROTACs.append(protac)
    except:
        fail_lnk = linker_df['smiles_dummy'].iat[idx] 
        fail_lnk_m = Chem.MolFromSmiles(fail_lnk)
        failed_linker.append(fail_lnk_m)
    idx += 1

print('Generated PROTACs : ', len(PROTACs))
print('failed linkers : ', len(failed_linker))

# 出力
# Generated PROTACs :  1092
# failed linkers :  7

1099個のリンカーのうち、1092個でユニットを繋ぐことができ、7個で失敗しました。

4-6. できたPROTACを確認してみよう!

それではどんなPROTACができたのか確認してみましょう。

最初の3つを取り出してみます。ちょっと描画を頑張ります。*4

from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

view = rdMolDraw2D.MolDraw2DSVG(990,350, 330, 350)
view.DrawMolecules([rdMolDraw2D.PrepareMolForDrawing(m) for m in PROTACs[:3]])
view.FinishDrawing()
svg = view.GetDrawingText()
SVG(svg)

f:id:magattaca:20220327210309p:plain

(-CH2CH2O-)の長さの異なるリンカーや、アミド結合の数が異なるものが出来ているようです。

もともとのARV-771のWarheadとリンカーの結合はアミド結合でした。ですので、リンカーライブラリの中で結合末端が炭素原子の場合は、Warheadとの結合がケトンになっています。あまり好ましくはないかもしれませんが、お試しということで。。。

「リンカーで多様性が出ているのか?」、記述子を計算して眺めてみましょう。

分子量MolWt)、LogPMolLogP)、水素結合アクセプター数CalcNumLipinskiHBA)、水素結合ドナー数CalcNumLipinskiHBD)を計算してデータフレームに格納します。*5

from rdkit.Chem import Descriptors, rdMolDescriptors

df_desc = pd.DataFrame(columns=['MW', 'LogP', 'HBD', 'HBA'])

df_desc['MW'] = [Chem.Descriptors.MolWt(m) for m in PROTACs] 
df_desc['LogP'] = [Chem.Descriptors.MolLogP(m) for m in PROTACs] 
df_desc['HBA'] = [Chem.rdMolDescriptors.CalcNumLipinskiHBA(m) for m in PROTACs] 
df_desc['HBD'] = [Chem.rdMolDescriptors.CalcNumLipinskiHBD(m) for m in PROTACs] 

ヒストグラムをプロットします。

import seaborn as sns
from matplotlib import pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(10,10))

sns.histplot(df_desc['MW'], kde=True, binwidth=50, ax=axes[0,0]) 
sns.histplot(df_desc['LogP'], kde=True, binwidth=0.5, ax=axes[0,1]) 
sns.histplot(df_desc['HBA'], kde=True, binwidth=1, ax=axes[1,0]) 
sns.histplot(df_desc['HBD'], ax=axes[1,1]) 

f:id:magattaca:20220327210356p:plain

思っていたよりも記述子にバラツキがでました。分子量では900 ~ 1500程度、LogPでも4 ~ 12前後となっています。

なお、今回もとにした化合物、ARV-771では以下の値となります。

print('MW of ARV-771 : ', Chem.Descriptors.MolWt(ARV771))
print('LogP of ARV-771 : ', Chem.Descriptors.MolLogP(ARV771))
print('HBA of ARV-771 : ', Chem.rdMolDescriptors.CalcNumLipinskiHBA(ARV771))
print('HBD of ARV-771 : ', Chem.rdMolDescriptors.CalcNumLipinskiHBD(ARV771))

# 出力
# MW of ARV-771 :  986.6619999999997
# LogP of ARV-771 :  6.531180000000008
# HBA of ARV-771 :  16
# HBD of ARV-771 :  4

リンカーライブラリと組み合わせて作ったPROTACライブラリと比較すると、一番頻度が高い範囲(より少し小さいくらい?)に含まれていそうです。

4-7. 失敗したリンカーを確認

最後に、molzipがうまくいかなかったリンカーを確認してみましょう。

view = rdMolDraw2D.MolDraw2DSVG(990, 600, 330, 200)
view.DrawMolecules([rdMolDraw2D.PrepareMolForDrawing(m) for m in failed_linker])
view.FinishDrawing()
svg = view.GetDrawingText()
SVG(svg)

f:id:magattaca:20220327210503p:plain

結合末端が、最初の5つは「窒素原子との2重結合( N=[* : 1])」、残りの2つは「同じAtomMapNumber([* : 2])が2つ」となっています。同じ結合箇所に2回結合を作ろうとして失敗しているようです。

エラーメッセージも「Explicit valence for atom # 18 C, 5, is greater than permitted」とでていました。「炭素の価数が5になっちゃうよ!」って怒ってます。

ごめんなさい。

5. まとめ

以上、今回はユニットで構成される化合物の例としてPROTACを題材に、RDKitのmolzip関数を使ってフラグメントを連結する方法で遊んでみました。また、PROTACで使用されるリンカー構造の取得元としてPROTAC-DBを参照しました。

生成したPROTACライブラリの記述子を計算すると、予想以上にバラツキがでました。WarheadE3 ligandと比べるとLinkerは単純な構造に見えますが、ユニットを組み合わせてできる化合物の性質には少なからず影響を与えそうです。

今回の記事では結合形成部位の化学的性質については全く考慮しませんでした。従って、記事中で見たように元々はアミド結合だったものがケトンになっていたりとしました。結合のタイプも含めて化合物生成を行うためには、RDKitのReactionといったものを使った方が良さそうです*6

ですが、Reaction SMILES(やSMARTS)はSMILESよりハードルが高くてしんどいので、お手軽に構造生成できるmolzipはとても便利だなーと思いました。

出てくるコードは基本的に参考にさせていただいた記事のコピペです。すみません。。。題材を変えたということで許してください。

おしまい*7

*1:このあたりおそらくSMILESではなくSMARTSの領域になると思うのですが、RDKitではMolFromSmilesで変換できます。表記方法に詳しくないのでよくわからないです。

*2:ググって一番上にでてきました。これが一般的な話かは分かりません。すみません。

*3:PROTAC-DBでARVで検索して一番上に出てきたのがこれでした

*4:参考:化学の新しいカタチさんの記事「RDKitでの構造式描画を詳しく解説

*5:参考:化学の新しいカタチさんの記事「RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ

*6:参考:化学の新しいカタチさんの記事「RDKitで化学反応:ケモインフォマティクスにおける反応式の扱い方

*7:ところで、今回の記事の冒頭にBuchwald-Hartwigクロスカップリングを上げてみました。色々なアミン芳香族をつなげて素晴らしい反応ですよね!アミンの種類が多ければ多いほどバラエティに富んだ構造をつくれます。 要するにビルディングブロックの多様性が、ケミカルスペースの探索空間につながるそうです。

Grygorenko, O.O., Volochnyuk, D.M., Vashchenko, B.V.
Emerging Building Blocks for Medicinal Chemistry: Recent Synthetic Advances
Eur. J. Org. Chem. 2021 6478-6510

OpenMMをステップバイステップで 〜 Part 2 PDBFixerでタンパク質の前処理 〜

前回の記事でMD計算に向けたタンパク質構造の前処理について、OpenMMに付随するGUIツールOpenMM Setupを使った操作を試しました。GUI操作の流れを把握する上では便利ですが、処理を自動化することを考えるとコマンドプログラムで扱えた方が良さそうです。

OpenMM SetupではPDB構造の修正に主にPDBFixerが使われています。また、TeachOpenCADDではPDBFixerのPython APIをつかうことで前処理を行なっています*1

そこで今回はJupyter notebook上でPDBFixerを使って前回の操作をたどってみます。といってもTeachOpenCADDのT019 · Molecular dynamics simulationからコードを切り貼りさせていただくだけですけどね!

f:id:magattaca:20220213124031p:plain
フィクサー*2

1. 処理の流れの確認

まずは処理の全体を再確認します。対象とするタンパク質の構造はPDB ID: 3POZで、以下のような処理を行いました。

f:id:magattaca:20220213124153p:plain
前処理の流れ

以上の操作をPythonベースで再現します。

2. PDBFixerで前処理を再現しよう

2-1. 準備

今回使うツールはOpenMMPDBFixerです。PDBFixerもcondaでインストールできます。

conda install -c conda-forge pdbfixer

タンパク質構造の処理のメインはPDBFixerですが、分子力場の設定水ボックスの配置、またPDBファイルの出力といった操作でOpenMMを利用することになります。*3

2-2. 処理の全体

最初に処理全体のコードをまとめておきます。TeachOpenCADDT019から抜粋して先に見た流れに合うように構成しています。*4

# 処理したいPDBファイルの指定
pdb_file = 'Data/3poz.pdb'

# 必要なライブラリのインポート
import openmm as mm
import openmm.app as app
from openmm import unit
import pdbfixer

# PDBFixerでPDBファイルを読み込む
fixer = pdbfixer.PDBFixer(pdb_file)

# 分子力場の設定
forcefield = app.ForceField("amber14-all.xml", "amber14/tip3pfb.xml")

# 不要な構造の削除
fixer.removeHeterogens()

# 欠けている残基のチェック(欠損原子の確認のためにも必要)
fixer.findMissingResidues()

# タンパク質末端の欠けている残基を取り除く処理
chains = list(fixer.topology.chains())
keys = fixer.missingResidues.keys()
for key in list(keys):
    chain = chains[key[0]]
    if key[1] == 0 or key[1] == len(list(chain.residues())):
        del fixer.missingResidues[key]

# 非標準な残基が含まれているか確認、あれば標準的なものに置き換える
fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()

# 欠けている原子の確認、あれば追加する
fixer.findMissingAtoms()
fixer.addMissingAtoms()

# 水素原子の付与(pHを設定する)
ph = 7.0
fixer.addMissingHydrogens(ph)

# 水ボックスの追加(力場、paddingの厚み、イオン濃度(デフォルトはNaCl))
modeller = app.Modeller(fixer.topology, fixer.positions)
modeller.addSolvent(forcefield, padding=1.0 * unit.nanometers, ionicStrength=0.15 * unit.molar)

# 処理後の状態(トポロジー、原子の位置)をPDBファイルで出力
top = modeller.getTopology()
pos = modeller.getPositions()
app.PDBFile.writeFile(top, pos, open('Data/processed.pdb', 'w'))

上のコードを実行することで、作業ディレクトリのData/3poz.pdbファイルを読み込んで処理した後、同じDataディレクトリに処理後のPDBprocessed.pdb)が書き出されます。

f:id:magattaca:20220213124248p:plain
五穀米おにぎり

PyMolで見た感じ再現できてそうです!

なお、追加したイオン(Na+、Cl-)の位置毎回変わることはご了承ください。*5

今回の記事の要点は以上で、あとは一つずつ中身を見ていくだけです。

3. ひとつずつ処理を確認

では、順番に操作の内容を確認してみましょう!

PDBFixerはデスクトップアプリコマンドラインアプリとしても使えるようですが、Python APIとして使うのが最も機能を有効利用できるらしく、マニュアルでもおすすめされています。

3-1. PDBFixerオブジェクト

PDBFixerをつかうには、処理したいPDBファイルを渡してPDBFixerオブジェクトを生成します。このオブジェクトに付随するメソッドを順番に呼んでいくことで、PDBファイルに対する処理が実行されていきます。

オブジェクトを作るのはこんな感じ。

fixer = pdbfixer.PDBFixer(pdb_file)
print(type(fixer))
#  <class 'pdbfixer.pdbfixer.PDBFixer'>

PDBFixerオブジェクトはトポロジーtopolgy)と位置positions)の組み合わせという形で構造情報を保持しています。positionsはそのまま原子の位置です。topologyもう一段上の抽象的な構造情報(ファイルに含まれるchain、各chainの残基、各残基の原子 etc.)といったものを指します。

試しにtopologyを表示してみます。こんな感じ。

print(fixer.topology)
# <Topology; 2 chains, 423 residues, 2536 atoms, 2463 bonds>

構造全体を要約したような情報が出てきました。

同様にpositionsも確認できます。

print(fixer.positions[0:3])
# [Vec3(x=-0.032400000000000005, y=3.2759, z=-0.39270000000000005), Vec3(x=-0.0237, y=3.1343, z=-0.3461), Vec3(x=0.119, y=3.0783000000000005, z=-0.3514)] nm

多いので最初の3つだけ表示しています。リストのようなオブジェクトなのでスライスで切り出せます。

上記の通り(x, y, z)の3次元座標のベクトルで各原子の位置を表しているようです。参考までにPDBファイル(3poz.pdb)をテキストエディタで開いて一つ目の原子の情報を確認すると以下のようになっています(各フィールドに注釈を加えています)。

f:id:magattaca:20220213124528p:plain

デカルト座標(⑦)の値を確認すると、PDBFixerオブジェクトのpositionsと一桁ずれています。PDBファイルはオングストローム(Å)、PDBFixerはナノメートル(nm)という単位の違いを踏まえれば一致していそうです。

PDBFixerオブジェクトがPDBファイル元にした構成であることがわかりました。*6

ちょっとくどくなりましたが、topologypositionで構造を表すという内容はOpenMMでも共通(?)で、処理後のPDBファイル出力の際にも関係してくるのでふれておきました。

3-2. 分子力場の設定

つづいて分子力場の設定です。「タンパク質構造の修正」という観点では関係なさそうですが、最後に「明示的な水分子を配置して水ボックスを作成する」という操作を行う上で必要になります。

ここではPDBFixerではなく、OpenMMappに含まれるForceField() を使います。OpenMMでは力場の情報がXMLファイルで定義されており、それを指定して読み込むことで利用できます。*7

以下のように力場を設定します。

forcefield = app.ForceField("amber14-all.xml", "amber14/tip3pfb.xml")

引数の1つ目メインの力場で、2つ目水モデルを指定する力場です。

メインで使っているamber14-all.xmlAMBER14力場の複数のファイルをまとめて読み込めるものです。種々の分子タイプ(amber14/protein.ff14SB.xmlamber14/DNA.OL15.xmlamber14/RNA.OL3.xmlamber14/lipid17.xml)をひっくるめて対応できて便利です。

水分子のモデルについてはamber14ディレクトリの中のamber14/tip3pfb.xmlを指定します。amber14ディレクトリの中のファイルは水だけでなく相当するイオンについての設定も含んでいるそうなので、このように力場を設定しておくほうが良いそうです*8*9

3-3. 異質分子(heterogen)の除去

つづいてタンパク質以外の異質分子(heterogen)の除去です。「PDB ID : 3POZ」では硫酸イオン(SO4)や低分子リガンド(03P, TAK-285)といったものを除きつつ、水分子(HOH)は残すといった処理を行います。

以下の1行でおしまいです。

fixer.removeHeterogens()

PDBFixerのコード(pdbfixer.py)を参照するとわかりますが、デフォルトの引数ではkeepWater=Trueとなっているため、水分子(HOH)は残されることになります*10

3-4. 欠けている残基の確認と除去(追加)

続いて欠けている残基(missing residue)の確認とそれらを「除去する or 追加する」を選択する操作です。ここで見つかったmissing residueは後ほど欠けている原子を追加する操作で補完されることになります。

以下のメソッドで欠けている残基を確認できます。

fixer.findMissingResidues()

見つかった欠損残基は以下のようにmissingResiduesで確認できます。(見やすさのためpprintで表示しています)*11

import pprint

pprint.pprint(fixer.missingResidues)
'''
    {(0, 0): ['GLY', 'GLU', 'ALA', 'PRO', 'ASN'],
     (0, 33): ['GLU', 'GLY', 'GLU', 'LYS'],
     (0, 43): ['ARG', 'GLU', 'ALA', 'THR', 'SER', 'PRO', 'LYS'],
     (0, 156): ['GLU', 'TYR', 'HIS', 'ALA', 'GLU', 'GLY', 'GLY'],
     (0, 285): ['GLU', 'GLU', 'ASP', 'MET', 'ASP', 'ASP'],
     (0, 293): ['ILE', 'PRO', 'GLN', 'GLN', 'GLY']}
'''

print('found ', len(fixer.missingResidues), 'missing regions')
#     found  6 missing regions

辞書型に格納されています。欠損残基として6つの領域が見つかりました。前回OpenMM Setupで確認した結果を再掲しておきます。

f:id:magattaca:20220213124921p:plain
missing residues

同じ結果が見つかっていそうです。辞書のKeyはタプル(chain.index, indexInChain)となっています。残基の位置は「鎖の中の位置(indexInChain)」で単純な残基番号(UniProtなどで確認できるもの)とは異なることにご注意ください。

見つかった欠損残基について、最終的な計算対象のモデル構造に含めたくない部分を取り除きます。タンパク質の両末端(N末端とC末端)を取り除き、鎖の中間部分(実験で見えなかったループ構造 etc.)を残したい場合は以下のようにします。

chains = list(fixer.topology.chains())
keys = fixer.missingResidues.keys()
for key in list(keys):
    chain = chains[key[0]]
    if key[1] == 0 or key[1] == len(list(chain.residues())):
        del fixer.missingResidues[key]

先に見た通りmissing residueの辞書のKeyは(chain.index, indexInChain)でした。後者IndexInChainについて0key[1] == 0)であれば欠損残基が鎖のN末端であることを意味します。逆に、鎖の残基数len(list(chain.residues()))と一致した場合は鎖のC末端が欠損していることになります。これらに該当する箇所が見つかったら「delで削除する」というのが上記です。

この操作を行ったのちmissingResiduesは以下のように変化しています。

pprint.pprint(fixer.missingResidues)

'''
    {(0, 33): ['GLU', 'GLY', 'GLU', 'LYS'],
     (0, 43): ['ARG', 'GLU', 'ALA', 'THR', 'SER', 'PRO', 'LYS'],
     (0, 156): ['GLU', 'TYR', 'HIS', 'ALA', 'GLU', 'GLY', 'GLY'],
     (0, 285): ['GLU', 'GLU', 'ASP', 'MET', 'ASP', 'ASP']}
'''

print('found ', len(fixer.missingResidues), 'missing regions')
# found  4 missing regions

N末端とC末端の2つが消えて、辞書の中身は4つの領域となりました。

取り除かなかった欠損残基(この段階でmissingResiduesに残っているもの)は、次々項(3-6.)でaddMissingAtomsを行う際にPDBFixerが補完してくれるようです。ミッシングループの構築に相当する操作ですね。

なお、今回は行いませんでしたが、もし「全てのmissing residueを取り除きたい」場合は以下のように空の字辞書で置換してあげればOKです。

fixer.missingResidues = {}

3-5. 非標準な残基の確認と置換

つづいて非標準な残基を確認し、置換する操作です。これは前回の記事で確認したOpenMM SetupGUI)の処理にはありませんでした。

今回対象としている 「3POZ.pdb」にはそもそも「非標準な残基」が含まれていないようなので、そのためかもしれません。一応該当するコードは以下です。

fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()

NonstandarResiduesfindしてreplaceする。PDBFixerのコードの文法がだいたいわかってきました。

なお、非標準なアミノ酸結晶学上の目的で加えられたりするようです(結晶化しやすくしたりする?)*12。タンパク質結晶学は専門でないのでよくわからないです。

3-6. 欠けている原子の確認と追加

つづいて欠けている原子missing atoms)の確認追加を行います。前々項(3-5.)で行った欠損残基missingResidues)の確認を踏まえた操作です。

まずは確認です。構造の中で欠けている重原子(heavy atom)を見つけます。

fixer.findMissingAtoms()

見つかった欠損原子は① missingAtomsと② missingTerminalsにわけて格納されます。

print(fixer.missingAtoms)
#  {}

print(fixer.missingTerminals)
#  {<Residue 292 (LEU) of chain 0>: ['OXT']}

missingAtomsで、② missingTerminalsには「292番 LeuのOXT」が格納されていました。

比較のため、前回OpenMM Setupで確認した結果を再掲しておきます。

f:id:magattaca:20220213125358p:plain

同じ重原子が見つかっていますね!

では、欠けている要素を追加します。addMissingAtomsメソッドにより「① missingAtoms、②missingTerminals、③missingResidues」の3つのフィールドに格納された欠損要素が追加されます。((メソッドはAtomsとなっていますがResiduesaddされるのがちょっとわかりにくいです))

fixer.addMissingAtoms()

引数にランダムシードを指定できることからわかるように、原子の追加の仕方(座標)にはランダムさがあります。

3-7. 水素原子の付与

さて、重原子について欠損の確認と追加が終わったので、今度は水素原子を付与します。この際pHを指定しますが、デフォルトの残基のpKaに従って水素原子の付与の仕方が決まるようになっています。

ph = 7.0
fixer.addMissingHydrogens(ph)

上では明示的に「pH = 7.0」の設定をしていますが、デフォルトが7.0なので空でも同じです。

3-8. 水ボックスの追加

最後にタンパク質構造を取り囲むように水ボックスを追加します。

OpenMMのappに含まれるModellerオブジェクトに変換した後で、addSolvent()で水分子を追加します。

# 水ボックスの追加(力場、paddingの厚み、イオン濃度(デフォルトはNaCl))
modeller = app.Modeller(fixer.topology, fixer.positions)
modeller.addSolvent(forcefield, padding=1.0 * unit.nanometers, ionicStrength=0.15 * unit.molar)

先に見た通りPDBFixerオブジェクトはtopologypositionsで構造情報を保持しています。この2つを取り出してopenmm.app.Modeller()に渡すことでModellerオブジェクトに変換できます。

addSolvent()の引数は「力場(forcefield)、水の層の厚み(padding)、イオン濃度(ionicStrength)」などです。*13順番に見ていきますが、先にOpenMM Setupの該当箇所を再掲しておきます。

f:id:magattaca:20220213125441p:plain

3-8-1. 分子力場

分子力場は先に「3-2. 分子力場の設定」で言及した通りです。addSolvent()ではデフォルトの水モデルが「TIP3Pmodel='tip3p')」となっています。力場の設定との互換性に注意したほうが良さそうです。

3-8-2. パディング

paddingではタンパク質を取り囲む水の層の厚みを指定します。

addSolvent()は基本的に「直方体の箱(rectangular box)に水を詰めていく」操作です。箱のサイズの指定方法には「① (x, y, z)ベクトルで各軸方向の長さを指定(boxVectors)」するか、「② タンパク質を中心としてその周り水の層厚さを指定して配置(padding)」するかの2通りあります。

「② padding」で指定する場合、タンパク質(溶質, solute)のxyz軸方向のうち一番大きい次元を求めて、その次元方向の両端にパディング距離分の厚みを加えた立方体の箱となります。つまり、1辺の長さ「最大次元 + パディング距離の2倍」の立方体です。

こういうイメージでしょうか?

f:id:magattaca:20220213125503p:plain

3-8-3. イオン濃度と中和

ionicStrengthで指定した濃度となるようにイオンが配置されます(単位 = モル濃度, molar)。

デフォルトではイオンはNaClです。「positiveIon='Na+'」「negativeIon='Cl-'」といったように引数でイオンの種類を設定できるので必要に応じて変更します。

選択できるイオンはそれぞれ「‘Cs+’, ‘K+’, ‘Li+’, ‘Na+’, ‘Rb+’」「 ‘Cl-‘, ‘Br-‘, ‘F-‘, ‘I-‘」です。OpenMM Setupで見たものと同じですね。

ところでイオンの追加には系の中和も関係してきます。デフォルトではneutralize=Trueとなっているので、溶質が電荷を帯びている場合には、その電荷を打ち消してボックス全体の電荷の合計がゼロとなるようにイオンが追加されます。

なお、水ボックスへの配置の手順は以下の通りです。

f:id:magattaca:20220213125523p:plain

上記の通り「ランダムに選んだ水分子をイオンと交換する」という配置の仕方なので、イオンの位置はランダムに変わります。

3-9. PDBファイルの出力

全ての処理が終わったので、最後に処理後の構造をPDBファイルで出力します。OpenMMappにあるPDBFile.writeFile()を使えばOKです。

構造情報としてtopologypositions、それから出力先のファイル名を引数にわたしてやります。

# トポロジーと位置情報をそれぞれ取得
top = modeller.getTopology()
pos = modeller.getPositions()

# 出力先ファイル名とともに引数に渡して書き出す。
app.PDBFile.writeFile(top, pos, open('Data/processed.pdb', 'w'))

これでprocessed.pdbというファイルがDataディレクトリの中にできました。

おしまい!

3-10. (おまけ)+αの情報

PDBFixerの使い方は以上ですが、上記で触れなかったOpenMMの使用法に関する情報を追加します。

3-10-1. openmmのunitライブラリ

水ボックスを設定する際にナノメートルモル濃度といった単位(unit)を利用しました。

ここでOpenMMのunitを利用しています。unitを使うと、数値とともにunit.nanometersunit.molarを記述するだけで、通常の単位を記述するのと同じ書き方で単位込みの情報として扱うことができます。

untiはシミュレーションを設定する際に時間(フェムト秒、ピコ秒、ナノ秒)を記述する際にも使えるので、OpenMMを利用する上でおさえておいたほうが良さそうです。

3-10-2. PDBFixerのaddSolvent

今回はTeachOpenCADDの記述に合わせてModellerオブジェクトに変更したあとでaddSolvent()を実行しましたが、PDBFixerにも同じメソッドは組み込まれているので、オブジェクトの変更が無駄に見えるかもしれません。

PDBFixerのコードを見ると、実際の操作はOpenMMを利用していて、「①Modellerオブジェクトに一度変換、②addSolvent、③PDBFixerオブジェクトに戻す」ということが行われています。

ここでは、水ボックス作成後、そのままOpenMMでPDBファイルを出力して終わりですので、PDBFixerオブジェクトに戻す操作がない分、記事のように記載したほうが無駄が少ないかもしれません。

4. 終わりに

以上、今回の記事では前回OpenMM Setupで行ったタンパク質構造の前処理操作を、PDBFixerPython APIで実行する方法を確認しました。

GUIで操作の流れを確認した後だと、だいたいやりたい処理の流れがわかるのがいいですね!

ひとつずつ処理を追ったせいで無駄に長くなってしまいました。毎度くどくてまとまりがなくてすみません。次回は、前処理済みの構造を使ったシミュレーションの設定方法を見ていきたいと思います。

今回も色々と間違いが多そうです。ご指摘いただければ幸いです。

ではでは!

*1:OpenMM Setupを利用していません

*2:イラストはいらすとやさんより

*3:PDBFixerは内部でOpenMMを利用しているので、OpenMMをインストールしていないと使えません。念のため・・・

*4:TeachOpenCADDはCC-BY-4.0ライセンスで公開してくださっています。感謝!

*5:後で確認しますが、設定した濃度(ionStrength)にあうをランダムに配置しています

*6:尚、3poz.pdbファイルには異方性温度因子ANISOU)の情報が各原子について記述されています。話が難しくなるので省略しています。

*7:力場についての詳細や利用できる力場の種類はドキュメントOpenMM User Guide : 3.6.2. Force Fieldsをご参照ください

*8:OpenMM User GuitdeのTIPの記載より

*9:尚、水分子の力場をTIP3P-FBとしていますが、このFBはForceBalanceの略のようです。元の文献はこちらです。(Lee-Ping Wang, Todd J. Martinez, and Vijay S. Pande Building Force Fields: An Automatic, Systematic, and Reproducible Approach J. Phys. Chem. Lett. 2014(5)1885DOI;10.1021/jz500737m

*10:他に残される残基は、通常の残基とUNK残基などです

*11:参考: note.nkmk.meさんの記事「Pythonのpprintの使い方(リストや辞書を整形して出力)

*12:PDBFixerのマニュアルの記述より

*13:引数の詳細はマニュアル addSolventを参照してください。