Zenn
🌟

脂質二重膜＋たんぱく質のMD

2023/05/09に公開
#
物理シミュレーション
#
gromacs
#
膜たんぱく質
tech

この記事は、Amberのチュートリアルのうち、Building and Equilibrating a Membrane System with PACKMOL-Memgen
でたんぱく質を埋め込んだ脂質二重膜を作り、脂質二重膜（たんぱく質なし）のチュートリアル「Simulating a Lipid Bilayer Tutorial」　を応用してデータを取るMDを実際に流してみたものになります。

ただ、ウチは一般のご家庭にあるパソコンでGROMACSしか使えません。なので後半は実行できるように都合よくパラメータ等を読み替えています。

お断り

たんぱく質とかやったことないけどMDシミュレーションなら多少知ってる程度の知識で書いています。
この記事の方法で得られる結果の保証はできません。

必要なもの

ソフトウェア

検証したシステム

  • WSL2（Ubuntu 22.04）
  • AmberTools 23（そのバージョンのアーカイブを展開するとamber22というフォルダができるやつです）
  • GROMACS 2022.4 （ビルド自体は GTX1060使っているときのもの）

ハードウェア

  • CUDAのコードが実行できるGPU
  • そこそこ早いCPU

ウチはi7-12700とRTX 3060を使いました。

モデル作成

この章はBuilding and Equilibrating a Membrane System with PACKMOL-Memgen　にあるコマンドをそのまま抜粋しただけになります。

ただこのページ、

  • 脂質二重膜のみ
  • たんぱく質を埋め込む
  • 溶質も加える

というそれぞれの段階でコマンドを紹介していますので、斜め読みだと不必要にコマンドを打つことになるかもしれません。

修正

AmberToolsのpackmol-memgen というものを使うのですが、pythonやnumpyのバージョンによってはエラーが出るかもしれません。
ウチは amber22/lib/python3.10/site-packages/packmol_memgen/lib/pdbremix にある v3numpy.py で、 np.float を np.float64 に書き換えてやる必要がありました。

準備

まず必要なファイルのダウンロード。

wget https://files.rcsb.org/download/1BL8.pdb
wget https://ambermd.org/tutorials/advanced/tutorial38/packing/ligand/TEA.gout
wget https://ambermd.org/tutorials/advanced/tutorial38/packing/ligand/leap_TEA.in

リガンド

こういう分子 TEA が周りにちょっとあるようです。これを脂質二重膜の周りの水層に少し溶かしてみます。
ＭＤシミュレーションで使えるようにするためには、

  • 何らかの量子化学計算で原子ごとの電荷を計算
  • 分子の結合情報からMDで使うパラメータを割り当て
  • AmberToolsで扱えるようなフォーマットに

という作業が必要です。最初の量子化学計算はチュートリアルのほうでやった結果ファイルを配布してくれているので、それを使います（準備で wget した TEA.gout です）。
以下のコマンドで処理できます。

antechamber -i TEA.gout -fi gout -o TEA_resp.mol2 -fo mol2 -c resp -rn TEA -pf y -at gaff2
parmchk2 -i TEA_resp.mol2 -f mol2 -s 2 -o TEA.frcmod
tleap -f leap_TEA.in

これは小分子 TEA について、

  • Gaussianによる量子化学計算した結果ファイル TEA.gout から RESP電荷を抽出、 TEA_resp.mol2に格納。
  • parmchk2 では TEA_resp.mol2 をもとに標準的な力場パラメータで対応できるかチェック、 無ければ推測して TEA.frcmod に追加。
  • tleap でMDに使える用のトポロジーファイルなどを作成。

これも脂質二重膜に組み込みます。

配置

packmol-memgen 一発です。

追記：残基名の編集

たんぱく質の周りに脂質を並べるだけなら一発なのですが、ちょっとMDを見据えると下処理が必要になります。
packmol-memgen 内部では tleap コマンドでパラメータ付けを行うのですが、どうもその時残基名が 1BL8.pdb のままでは不都合があるようです。
特にpH調整の結果イオン化したり中性化したりしたらそれ用に変えてやる必要があるようです。

ということで、まずパラメータ無しで構造PDBファイルのみを生成します。

packmol-memgen --pdb 1BL8.pdb --lipids DOPE:DOPG --ratio 3:1 --keepligs --solute TEA.pdb --solute_con 4 --solute_prot_dist 10

これで bilayer_1BL8_lipid.pdb ができます。で、残基名変更は完全に対象により違う操作になります。 今回のケースではGLU を GLH に置換します。

for i in " 49" 146 243 340; do sed -i "/GLU . $i/s/GLU/GLH/g" bilayer_1BL8_lipid.pdb; done

なお、手元で試したときは忘れてたかもしれない。

あとは追記前に合流。

パラメータ付け

packmol-memgen --pdb 1BL8.pdb --lipids DOPE:DOPG --ratio 3:1 --keepligs --solute TEA.pdb --solute_con 4 --solute_prot_dist 10 --ligand_param TEA.frcmod:TEA.lib --gaff2 --parametrize

先のコマンドで bilayer_1BL8_lipid.pdb ができていたら、 packmol は終わったものとしてスキップ、続きの tleap でのパラメータ付けを行ってくれます。

軽くオプションの意味を解説します。

Option 引数 解説
--pdb 1BL8.pdb
--lipids DOPE:DOPG 2種類の脂質を膜に使う
--ratio 3:1 DOPE:DOPG の存在比
--keepligs --pdb で渡された構造に付随するリガンドを消さない
--solute TEA.pdb 追加で一緒に溶かす物質
--solute_con 4 溶質の濃度
--solute_prot_dist 10 たんぱく質からこの距離内の円柱の中に溶質分子を配置
--ligand_param TEA.frcmod:TEA.lib 溶質小分子のパラメータファイル
--gaff2 gaff2 でパラメータを割り当て
--parametrize MD用パラメータファイルを生成する

およそ10分ぐらいかかりました。

これをやった結果が bilayer_1BL8_lipid のbasenameに拡張子 .crd .pdb .top で3つのファイルができているかと思います。

この時点でのpdbはこんな感じ

変換

Amber用のファイルをGROMACSに変換します。以下のPythonコードを同じディレクトリで実行しましょう。

import parmed as pmd
parm = pmd.load_file('bilayer_1BL8_lipid.top', 'bilayer_1BL8_lipid.crd')
parm.save('system.top', format='gromacs')
parm.save('bilayer_1BL8_lipid.gro')

すると system.top と bilayer_1BL8_lipid.gro の2つのファイルができます。

MD 実行

チュートリアルとの差異

ここからは「Simulating a Lipid Bilayer Tutorial」　の内容を追ったものになります。
ただAmberのチュートリアルをGROMACSでやる都合上、あと個人の好み、たんぱく質が増えるという違いからくる試行錯誤の結果で少々異なる部分があります。

  • SHAKEではなくLINCS
  • BarostatはGROMACSのWarningが出たのでParrinello-Rahman
  • チュートリアルではHoldというステップがありますがこれはAmberではGPU版のほうで相互作用原子リストを更新しないため、 skinnb Error というのが生じるそうで、その対策だそうです。理解した限りではGROMACSでは nstlist で適切に設定すれば解決する問題かなと思いますので、Holdステップはスキップし、長めの平衡化MDとしました。

下処理

個人の好みですが、シミュレーションの系の箱のサイズで X と Y が同じになるようにして、安心して Semiisotropic な熱浴を適用できるようにしてもいいかなと思います。
gro ファイルの一番下は

11943WAT     H263465   3.411   6.494   8.458
   9.29350   9.27380   9.84606

となっていると思いますが、この一番下の行の3つの数字で最初の2つの数字を大きいほうに合わせます。
このケースだと 9.29350 9.29350 9.84606 という感じですね。
物理的には隙間ができることに相当しますが、そもそも発生させた構造の時点で空いているので問題ないかなと思います。

なんか最初から X = Y となるようなオプションとか packmol-memgen にないかな

あと、MD実行時に初めて出現するエラーを見ないと分からない問題でうまくいかないことがあります。
これはファイルを見ただけではわかりませんので、いったん最初のエネルギー最小化をやってみましょう。
ここでエラーが起こらず（力が0でない原子がある、程度だったら目をつむります） 最小化の結果の構造である min.gro ができればOKですが、もしエラーが出た場合、

  • エラー・Warningが出ない構造になるまで先の packmol-memgen をリセマラ
  • 適宜 gro か top ファイルを修正

のどちらかを行います。以下は考えられるエラーの原因2通りで、これに合わせて修正しましょう。なお両方起こる可能性もあります。

電荷が0でない

gmx grompp でWarningの形で出ます。 トポロジーファイルを見ていると残基ごとに qtot を計算しています。
普通は 1.00000 とか 0.000000 といった整数なのですが、これがたまに 1.0000003 みたいに誤差がついてて、これが合計したときに消えない時があります。
こういう時は少々電荷が変わっても良さそうな原子（対称的な原子がないやつとかで変化量が0.1% 未満になりそうなもの）を書き換えて帳尻を合わせてみます。

原子の位置が近すぎる

gmx mdrun で

# Steepest Descents converged to machine precision in 11 steps,
# but did not reach the requested Fmax < 10.
# Potential Energy  =  1.0547430e+17
# Maximum force     =            inf on atom 16811
# Norm of force     =            inf

みたいな形で力が発散したというエラーで止まります。
このときは該当原子をgroファイル中で 0.01 程度動かすとエラーが解消されました。

エネルギー最小化

まずインプットファイル min.mdp を作ります。

min.mdp
integrator   = steep
nsteps       = 20000
emstep = 0.0002
dt           = 0.0002
constraints              = h-bonds
constraint-algorithm     = Lincs
continuation             = no
Lincs-iter = 1
Lincs-order = 4
nstlog       = 1
nstlist      = 1
ns_type      = grid
pbc          = xyz
rlist        = 1.0
coulombtype  = PME
rcoulomb     = 1.0
rvdw         = 1.0
cutoff-scheme = Verlet

これをコンパイルして実行します。

gmx grompp -c bilayer_1BL8_lipid.gro -p system.top -f min.mdp -o min.tpr
gmx mdrun -deffnm min

ちなみに試行錯誤の過程でLINCSを付けてない最小化 → つけた最小化という2段階のMDをやっているのですが、見る限り付けてない最小化は要らないかなと思いここでは省略しました。

平衡化

以降、Productionも含めすべてのMDはNPTです。

100 Kに

最初の系は絶対零度ですので、いきなり300Kにすると壊れるかもしれません。まず100Kに暖めます。

npt1.mdp
integrator   = sd
nsteps       = 5000
nstxout      = 10
dt           = 0.002
constraints              = h-bonds
constraint-algorithm     = Lincs
continuation             = no
Lincs-iter = 1
Lincs-order = 4
cutoff-scheme = Verlet
nstlog       = 1
nstlist      = 10
ns_type      = grid
pbc          = xyz
rlist        = 1.0
coulombtype  = PME
rcoulomb     = 1.0
rvdw         = 1.0
tc-grps      = System
tau-t        = 1.0
ref-t        = 100
pcoupl       = Parrinello-Rahman
tau-p        = 1.0
ref-p        = 1.0 1.0
pcoupltype   = semiisotropic
compressibility = 4.5e-5 4.5e-5

pcoupltype を semiisotropic とすることでXY軸とZ軸で違う取り扱いとしました。試行錯誤のときはisotropic でもまあ問題なさそうでしたが、一応。

あと、熱浴の設定ですが、チュートリアルはAmberのインプットを見ると

gamma_ln=1.0,   ! Collision frequency for Langevin thermostat

と \gamma = 1.0 ps のランジュバン熱浴です。GROMACSではランジュバン熱浴を指定するとき、実は tcoupl ではなくintegrator = sd で指定します。また tau-t は \gamma の逆数(1/ps) を指定します。
これで上で示したようなインプットパラメータになります。 integrator が md とか md-vv でないのはこのためです。まあ普通に tcoupl = V-rescale とかでも十分だと思います。

これをコンパイルして実行します。

gmx grompp -c min.gro -p system.top -f npt1.mdp -o npt1.tpr
gmx mdrun -deffnm npt1

これで 100K に温まった構造 npt1.gro が得られるかと思います。

300 Kに

１からインプットファイルを作るのはアレなので、書き換える形で作成します。

sed 's/nsteps *= *5000/nsteps = 50000/;s/nstxout *= *10/nstxout = 500/; s/ref-t *= *100/ref-t = 300/;s/tau-p *= 1.0/tau-p = 2.0/' npt1.mdp > npt2.mdp

変更したところを抽出すると以下のようになります。

nsteps = 50000
nstxout = 500
ref-t = 300
tau-p = 2.0

nstxout を変えているのは、出力された trr ファイル（軌跡が格納されます）が必要以上に大きくならないようにするためです。
これでMDを実行します。

gmx grompp -c npt1.gro -p system.top -f npt2.mdp -o npt2.tpr
gmx mdrun -deffnm npt2

最後の仕上げ

最後に長めに 300 K, 1 atm の平衡化をしておきます。

sed 's/nstlog *= *[0-9]*/nstlog = 10000/;s/nsteps = *[0-9]*/nsteps = 2500000/;s/nstxout = *[0-9]*/nstxout = 25000/;s/tau-p = 2.0/tau-p = 1.0/' npt2.mdp > npt3.mdp

変わったところの確認。

nsteps = 2500000
nstxout = 25000
nstlog = 10000
tau-p = 1.0

シミュレーション時間が長いと、情報の少ないLogファイルもアホみたいにでかくなるので、 nstlog も変更します。

gmx grompp -c npt2.gro -p system.top -f npt3.mdp -o npt3.tpr
gmx mdrun -deffnm npt3

ProductionMD

インプットファイルは最後の npt3.mdp を使いまわします。

sed 's/nsteps = [0-9]*/nsteps = 2500000/;s/nstxout = [0-9]*/nstxout = 2500/;' npt3.mdp > prod.mdp

変わったところを確認。

nstxout = 2500

実は軌跡を出力する頻度 nstxout が10倍になっただけで nsteps が同じでした。

gmx grompp -c npt3.gro -p system.top -f prod.mdp -o prod.tpr
gmx mdrun -deffnm prod

これでできた prod.trr や prod.edr を使って解析していきます。

結果確認

動画

実際に原子たちがどう動いてみるかを確認します。
まず軌跡データですが、そのままだと周期境界条件を気にせず座標を時間発展しているので、そのままVMDで出力するとどこかに拡散していってしまいます。
実際の力の計算のときは周期境界の長さで割り算したあまりを使うので問題はないですが、見た目の問題です。

そこで周期境界に合わせた出力を行います。今回はたんぱく質の位置などにこだわりがないので、ほぼ標準的なコマンドで大丈夫です。

gmx trjconv -f prod.trr -o prod_nop.trr -s prod.tpr -pbc whole

どのグループを返還するか聞かれるので、 0 を選択します。

その後 npt3.gro をVMDに読み込み、その系に対しLoad data into Molecules... から prod_nop.trr を読み込みます。
で、 Representation を少々工夫するとこのような動画が見れるかと思います。

https://www.youtube.com/watch?v=yjVen7y8rpA

なお63465原子の1001ステップ分のデータで 728 MB ありましたので、処理が終わり次第削除しています。

電子密度

Z軸に沿った電子の密度を見てみます。
といっても量子化学計算をやるわけではなく、原子番号 Z の周りには一定の半径の中に Z 個の電子があるとして単純に足していったものになります。
gmx density コマンドにその機能があります。
ただGROMACSにはユーザーが好き勝手に付けたラベルと恣意的に決められる質量しか与えられていないので、各粒子が持つ電子の数を与える必要があります。
それを渡すには一定のフォーマットのファイル electrons.dat を作成します。

electrons.dat
167
C = 6
H = 1
K = 19
（略）

最初の 167 は何種類の原子ラベルがあるかの宣言、あとは各ラベルが何個の電子があるかを示します。原子ラベルはgroファイルのフォーマットだと 10～15文字目にあるもので、右詰めスペース埋めとなっています。これから重複を削除して列挙、それぞれについて最初の1文字を元素だと思って原子番号＝電子数を割り当て、 = 原子番号 と追加します。
さすがにこの量は手作業ではしんどいので、Bashコマンドやawkを活用します。

cut -c 10-15 npt3.gro | sort -u | awk '
BEGIN{
  nrs["H"]=1;
  nrs["C"]=6;
  nrs["N"]=7;
  nrs["O"]=8;
  nrs["P"]=15;
  nrs["S"]=16;
  nrs["Cl"]=17;
  nrs["K"]=19
}{
  sub(/^ */, "");
  elem=substr($0,1,1);
  print $0 " = " nrs[elem]
}' > electrons.dat

これで作る electrons.dat には最初と最後の行に１，２行ほど余計なものがあるのでそれを削除、行数を数えて1行目にその数字を入れます。

あ、記事にまとめている今気づいたけどバグがありますね、 元素名を1文字でしか見ていないので Cl が6になっちゃう。
まあ7個だけだし今回はほとんど結果に影響しないと思いますので放置します。
ただ2文字の元素が中心になってくると組み直さないといけないですね。

 gmx density -s prod.tpr -f prod.trr -dens electron -ei electrons.dat

これで系全体に対応する 0 を選択します。
問題なければ density.xvg が出力されますので、 xmgrace などでプロットします。

チュートリアルにある結果に比べて中央の構造が崩れていますが、これはたんぱく質が入ったことによるものかなと思っています。
たんぱく質 1BL8 は円錐台のような形をしていてZ軸が大きいほうがより太くなっているように見えました。それに合わせてたんぱく質が大きいほうがより崩れているかなと。

参考

packmol

packmol は分子のファンデルワールス半径などを考慮して、いい感じに分子を並べていくプログラムです。独特の文法を持っていて、たとえば

structure PROT0.pdb
  number 1
  fixed 0. 0. 0 0. 0. 0.
  radius 1.5
end structure

は 「PROT0.pdb の中身を1個、原点においてください」 という意味で、

structure DOPE.pdb
  nloop 20
  number 75
  inside box -41.84 -41.84 -23.0 41.84 41.84 0.0
  atoms 1 11
    below plane 0. 0. 1. -18.0
  end atoms
  atoms 79 126
    over plane 0. 0. 1. -4.0
  end atoms
end structure

は 「 DOPE.pdb は試行20回で、-41.84 < x, y < 41.84、-23 < z < 0 の範囲に75個おいてね。あと、1番目の原子と11番目の原子は 0*x + 0*y + z = -18.0 の平面の下に、79番目と126番目の原子は z = -4 の平面より上においてね」
という指示になります。

頑張れば球形ミセルなどいろんなものができるかなと思います。

初めは独立したプログラムだったのですが、ある時からambertoolsに組み込まれました。

VMD

今回の可視化はこういう感じのものです。

とりあえず bilayer_1BL8.pdb をVMDで読み込み、Graphics > Representations... で

という感じ（Materialは最初以外はOpaque）にすると、

という感じに見えます。

ちょっと苦労した点として、 K+ とか Cl- みたいに記号が入った名前を指定するときはクォーテーションマークで囲う必要があることでした。

パフォーマンス

最後のデータ取得MDでかかった時間です。
初めのしばらくはYouTubeとか見ながらやってたので、MDに集中させれば100 ns/day いけるかもしれません。

               Core t (s)   Wall t (s)        (%)
       Time:    46070.183     4607.022     1000.0
                         1h16:47
                 (ns/day)    (hour/ns)
Performance:       93.770        0.256

続き

https://zenn.dev/meguru_compchem/articles/28e54e3f6e5e45

Discussion