脂質二重膜＋たんぱく質のMD
この記事は、Amberのチュートリアルのうち、Building and Equilibrating a Membrane System with PACKMOL-Memgen
でたんぱく質を埋め込んだ脂質二重膜を作り、脂質二重膜（たんぱく質なし）のチュートリアル「Simulating a Lipid Bilayer Tutorial」 を応用してデータを取るMDを実際に流してみたものになります。
ただ、ウチは一般のご家庭にあるパソコンでGROMACSしか使えません。なので後半は実行できるように都合よくパラメータ等を読み替えています。
お断り
たんぱく質とかやったことないけどMDシミュレーションなら多少知ってる程度の知識で書いています。
この記事の方法で得られる結果の保証はできません。
必要なもの
ソフトウェア
- Linux もしくは WSL2（GPUが使えたほうがいいです）
- AmberTools
- GROMACS
- ParmEd
検証したシステム
- 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版のほうで相互作用原子リストを更新しないため、
skinnbError というのが生じるそうで、その対策だそうです。理解した限りでは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 を作ります。
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に暖めます。
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
と tcoupl ではなくintegrator = sd で指定します。また tau-t は
これで上で示したようなインプットパラメータになります。 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 を少々工夫するとこのような動画が見れるかと思います。
なお63465原子の1001ステップ分のデータで 728 MB ありましたので、処理が終わり次第削除しています。
電子密度
Z軸に沿った電子の密度を見てみます。
といっても量子化学計算をやるわけではなく、原子番号
gmx density コマンドにその機能があります。
ただGROMACSにはユーザーが好き勝手に付けたラベルと恣意的に決められる質量しか与えられていないので、各粒子が持つ電子の数を与える必要があります。
それを渡すには一定のフォーマットのファイル 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
続き
Discussion