您的位置 首页 Protocol

蛋白质聚集的分子动力学模拟

蛋白质聚集的分子动力学模拟

( Molecular dynamics simulations of protein aggregation)

蛋白质紊乱和聚集在许多神经退行性疾病(如阿尔茨海默病和帕金森病)的发病机制中起着重要作用。这些疾病中聚集过程的最终产物是高度结构化的淀粉样纤维。因此,合理设计靶向淀粉样蛋白寡聚体形成的药物,就需要全面了解驱动蛋白质聚集的物理化学力。原子级的分子动力学 (MD) 模拟提供了该过程最高的时间和空间分辨率,可以捕捉淀粉样蛋白寡聚体形成过程中的关键步骤。

1.MD 模拟和基本分析

以 Aβ16-22 为例,其两端都带有封端基团,序列为 ACE-KLVFFAE-NME。由于蛋白质数据库 (PDB) 不包含该肽的结构,因此可以从 Aβ42 的 PDB 结构的残基 16-21 的坐标中检索以下模拟的起始结构,例如 PDB ID 1Z0Q。在 PyMOL 的蛋白质模式下使用 Builder 工具,ACE 和 NME 封端组可以分别添加到 N 和 C 端。

1.1 六肽模拟体系建立

1. 第一步是为 Aβ16-22 单体产生一个松弛的构象。对于 Aβ16-22,建议模拟长度为 1 μs 或更长,使用构象聚类确定最稳定的单体结构(本文省略这步)。

2. 使用 PACKMOL将六个单体随机放置到一个模拟盒子中。下面的示例脚本将六个 Aβ16-22 肽放置在一个大小为 ~10 nm × 10 nm × 10 nm 的模拟框中,它们之间的距离至少为 1.2 nm(或脚本中的 12 Å)。
packmol.inp
#Six monomers of abeta16-22 peptide 
#minimum distance between two monomers 
tolerance 12.0 
seed -1 
#The file type of input and output files is PDB 
filetype pdb 
#The name of the output file 
output abeta16-22_hexamer.pdb 
#add TER to after each monomer chain 
add_amber_ter 
#distance from the edges of box 
add_box_sides 1.0 
#path to input structure file 
#units for distance is measured in Angstrom 
#box size is 100 Å 
structure abeta16-22.pdb 
number 6 
inside box 0. 0. 0. 100. 100. 100. 
end structure
运行packmol

 packmol < packmol.inp

或者使用gmx insert-molecules命令

gmx insert-molecules -ci abeta16-22.pdb -nmol 6 -box 10 10 10 -o abeta16-22_hexamer.pdb

1.2 不同的模拟步骤创建目录

创建五个主要步骤的目录:拓扑文件构建、能量最小化、NVT 平衡、NPT 平衡和 MD 生产运行。
每一步都需要一个 .mdp 文件。mdp 文件类型扩展名代表分子动力学参数,这些文件包含使用 GROMACS 设置 MD 模拟的所有关键参数。

mkdir 1-topol 2-em 3-nvt 4-npt 5-md mdp

1.3 拓扑构建

拓扑文件包含有关将被模拟的分子类型和分子数量的信息。上一步的 .pdb 文件作为输入,除了拓扑文件之外,还会生成一个 .gro 文件,该文件与 .pdb 文件一样,也包含模拟系统的坐标。它们之间的主要区别在于它们的格式。

(1)从http://macerell.umaryland.edu/download.php?filename=CHARMM_ff_params_files/charmm36-mar2019.ff.tgz下载Charmm36m力场,复制到1-topol目录下。切换到该目录:

cd 1-topol/

(2)运行GROMACS pdb2gmx 命令处理输入的结构文件并创建扩展名为.top 的拓扑文件、扩展名为.itp 的拓扑包含文件和扩展名为.itp 的位置约束文件。

gmx pdb2gmx -f ../abeta16-22_hexamer.pdb -o protein.gro -p topol.top -ignh -ter <<EOF 

EOF

选项说明:
-f: 读取输入结构文件 abeta16-22_hexamer.pdb
-o and -p: 写入输出结构文件 protein.gro 和系统拓扑文件 topol.top
-ignh: 忽略输入文件中的氢原子,这是可取的,因为输入文件和力场中氢原子的命名约定不同。GROMACS 将使用所选力场的 H 原子名称添加新的氢原子。
-ter: N 端和 C 端的质子化状态可以-ter 标志交互选择
Option 1: choosing protein force field (charmm36-mar2019)
Option 1: choosing water force field (TIP3P)
Option 3: choosing N-terminus (None, as we use ACE capping)
Option 4: choosing C-terminus (None, as we use NME capping)
Options 3 and 4 are repeated for each peptide in the system, in this example six.

成功执行 GROMACS pdb2gmx 命令后,该目录还将包含六个拓扑包含文件和六个位置约束文件,每个肽一个:topol_Protein_chain_$chain.itp 和 posre_Protein_chain_$chain.itp,其中 $chain = A、B、C、D、E 或 F。后面的文件包含肽的所有非氢原子的位置限制条目,这是平衡 MD 步骤期间需要的。

(3) 接下来,创建一个模拟框。请注意,上面定义的框仅用于放置肽。和以前一样,选择了一个 10 × 10 × 10 nm3 的立方体盒子:

gmx editconf -f protein.gro -o box.gro -bt cubic -box 10 10 10

(4)在模拟盒里装入水分子

gmx solvate -cp box.gro -cs spc216.gro -o protein-solvated.gro -p topol.top
溶剂分子的添加反映在 topol.top 文件中,除了 6 个肽段外,该文件现在还包括 32,323 个水分子。
(5)为了使用周期性边界条件进行 MD 模拟,我们需要根据需要通过添加正 (Na+) 或负 (Cl-) 离子来中和系统的电荷。此外,我们可以为系统指定一个特定的离子浓度,通常约为 150 mM 以模拟生理条件:

gmx grompp -f ../mdp/ions.mdp -c protein-solvated.gro -p topol.top -o protein-ions.tpr 

echo 13 | gmx genion -s protein-ions.tpr -o protein-ions.gro -p topol -neutral -conc 0.15 

Explanation of flags and options: 
-neutral: 通过添加 Na+ 或 Cl- 离子将系统中和至零电荷 
-conc: 将浓度更改为 150 mM NaCl
echo 13: 用离子替换溶剂分子,其中选项 13 代表 SOL(溶剂分子)。
在第一个命令中,grompp 模块,读取坐标文件、系统拓扑文件、ions.mdp 文件并将它们处理成 GROMACS 二进制格式生成 .tpr 文件。该文件包含模拟的起始结构、单个肽和水分子的拓扑信息以及所有模拟参数。第二个命令读取二进制输入文件 protein-ions.tpr 并中和净电荷,添加了 90 个 Na+ 和 90 个 Cl- 离子,相应地更新 topol.top 文件。

图 1 六个 Aβ16-22 肽(显示为不同颜色),随机放置在水盒子中

上面我们准备了模拟的拓扑文件和坐标文件,接下来运行MD模拟。

1. 能量最小化

(1)对于能量最小化 (EM) 步骤,切换到目录 2-em:

cd ../2-em

(2)如 em.mdp 文件所示,我们采用最速下降法来最小化系统,直到最大力为达到 100 kJ/mol/nm 或 2,000 个最小化步骤。grompp 用于将结构、拓扑和仿真参数组合成二进制输入文件 protein-em.tpr,然后将其传递给 GROMACS mdrun 命令。

gmx grompp -f ../mdp/em.mdp -c ../1-topol/protein-ions.gro -p ../1-topol/topol.top -o protein-em.tpr

gmx mdrun -v -deffnm protein-em 

Explanation of flags: 
-v: 将 EM 步骤的进度打印到屏幕上 
-deffnm: 定义输入和输出文件

成功执行 mdrun 命令后,会生成以下文件:

protein-em.gro: 最终能量.最小化结构文件 
protein-em.edr: 二进制格式的能量文件 
protein-em.trr: 轨迹文件,包括所有坐标、速度、力和能量的二进制格式 
protein-em.log  : ASCII 格式的 EM 步骤的文本日志文件

2. NVT 平衡 

在 EM 步骤之后,执行两个平衡步骤。第一个平衡步骤在等温和等容条件下进行,称为 NVT 系综,因为 N(粒子数)、V(体积)和 T(温度)保持恒定。此外,在此步骤中,只有溶剂分子和离子在肽周围达到平衡,使它们达到所需的温度(在我们的例子中为 300 K),而肽原子的位置受到限制。为此,使用包含在拓扑构建期间生成的位置约束的 .itp 文件。类似于 EM 步骤,在 NVT 和以下 MD 步骤 grompp 和 mdrun 被调用:

cd ../3-nvt 

gmx grompp -f ../mdp/nvt.mdp -c ../1-topol/protein-em.gro -p ../1-topol/topol.top -o protein-nvt.tpr -r ../2-em/protein-em.gro

gmx mdrun -v -deffnm protein-nvt

通常,NVT 平衡步骤是 100 ps 长的 MD 模拟,足以在所需温度下平衡蛋白质或肽周围的水。以下是 nvt.mdp 文件中设置的重要参数的说明:
gen_seed = -1: takes as random number seed the process ID. 
gen_vel = yes: generates random initial velocities. For this, the random number seed is used. 
tcoupl = V-rescale: defines the thermostat. 
pcoupl = no: pressure coupling is not applied.
如果 grompp 根据其进程 ID 分配随机数种子,则每次重新运行 grompp 时,它将分配不同的种子,因为该 grompp 执行的相应进程 ID 也不同。这保证了种子总是随机的,因此每次重复模拟时,都会生成不同的随机数,从而导致不同的初始速度分布。当模拟重复多次以收集系统统计数据时,这一点很重要。温度耦合是使用速velocity rescaling热浴,这是一种改进的 Berendsen 弱耦合方法。成功执行 mdrun 命令后,将生成具有与 EM 步骤中生成的文件类型扩展名相同的文件。

3. NPT 平衡

在第二个平衡阶段,系统的压力和密度使用等温等压系综进行调整,也称为 NpT 系综,因为 N、p(压力)和 T 保持不变持续的。

cd ../4-npt 

gmx grompp -f ../mdp/npt.mdp -c ../1-topol/protein-nvt.gro -p ../1-topol/topol.top -o protein-npt.tpr -r ../3-nvt/protein-nvt.gro

gmx mdrun -v -deffnm protein-npt

与 NVT 平衡运行相比,npt.mdp 文件中的一项更改是增加了压力耦合部分,使用 Parrinello-Rahman压浴。其他值得注意的变化是:
continuation = yes: continuation of the simulation from the NVT EQ step. 
gen_vel = no: velocities will be read from the trajectory files generated from the NVT equilibration step and not newly initiated.
平衡阶段成功执行后,整个系统的温度和压力都将得到调整,以便我们可以继续执行生产运行。

4. MD 生产运行 

在生产运行中,去除蛋白质的位置限制;然而,所有键长都将使用LINCS 方法限制在其平衡值,使用 2 fs 的时间步长对运动方程进行积分。为了对构象空间进行充分采样,我们需要以微秒的数量级运行生产运行。

要执行本例中的MD 生产模拟,切换到相应目录

cd ../5-md 

gmx grompp -f ../mdp/md.mdp -c ../1-topol/protein-npt.gro -p ../1-topol/topol.top -o protein-md.tpr -r ../4-npt/protein-npt.gro

gmx mdrun -v -deffnm protein-md

在 MD 生产运行之后,将生成与前面步骤相比的另一种文件类型,一个 .xtc 轨迹文件。 
在成功完成 MD 生产运行后,可以开始对 MD 轨迹进行分析。protein_md.xtc 文件包含由生产模拟采样的系统的所有坐标,用于分析。各种 GROMACS 分析工具可以读取 .xtc 文件的二进制格式。但是,在下文中,我们将使用基于 Python 的 MDAnalysis 和 MDTraj 之类的工具,这些工具也可以处理 .xtc 文件。接下来主要分析寡聚体大小和聚集体中的接触。

1. Anaconda创建环境

conda create -n conda-python3.6 python=3.6

conda activate conda-python3.6

安装MDanalysis 和MDtraj

conda install -c conda-forge mdanalysis

conda install -c omnia mdtraj

2. 为分析创建一个新目录

cd ../ 

mkdir analysis 

cd analysis/

MD轨迹处理

(1)将轨迹文件 protein_md.xtc 和运行输入文件 protein_md.tpr 从生产 MD 目录复制到分析目录。

cp ../5-md/ protein_md.xtc ./

cp ../5-md/ protein_md.tpr ./

(2)对于分析,我们只需要蛋白质的坐标,而不需要溶剂和离子的坐标。使用 GROMACS trjconv 命令执行提取和重新保存:

gmx trjconv -s protein_md.tpr -f protein_md.xtc -o protein_only.xtc 
On the command prompt: select 
option “1” for centering “protein”, 
option “1” for output “protein” 
The output .xtc file will include only coordinates of the peptide chain.

(3)我们还需要以 .gro 或 .pdb 格式提取参考结构文件。在当前示例中,我们创建了一个 .pdb 文件。

gmx trjconv -s protein_md.tpr -f protein_md.xtc -o protein_only.pdb -dump 0 

On the command prompt: select
option “1” for output “protein”
Explanation: -dump 0: 储轨迹文件的第一帧。

(4)需要在分析之前解决的蛋白质聚集模拟的一个特殊问题是 MD 模拟过程中的 PBC,它可能导致蛋白质似乎被破坏。许多分析脚本无法处理此类破碎的蛋白质,这会导致分析中出现假象。在 VMD 中读入新创建protein_only.pdb 和 protein_only.xtc 文件并可视化轨迹。然后在“Extensions”选项卡上打开 Tk 控制台,并通过输入命令可视化蛋白质系统周围的 PBC 框:

pbc box

[atomselect top all] set chain 0

pbc join fragment -all

另存为:protein-nopbc.trr

或者使用trjconv处理PBC

gmx trjconv -s protein_only.pdb -f protein_only.xtc -pbc nojump -o protein_nopbc.trr

4. 计算低聚状态和残基间接触频率

conda activate conda-python3.6

python oligos-cmap.py protein_only.pdb protein_nopbc.trr 4

python plot-cmap.py protein_only.pdb protein_nopbc.trr contact-map.dat

python plot-oligostate.py protein_only.pdb protein_nopbc.trr oligo-highest-size.dat 100

分析的结果如图 2 所示。图 2A 中显示六个肽在约 80 ns 开始形成五聚体。对组成寡聚体的肽之间的残基-残基接触进行计数,形成图 2B 中的概率图。它显示肽更喜欢以反平行方向组装,这通过带相反电荷的 N 端 Lys16 和 C 端 Glu22 残基之间的静电相互作用稳定。此外,还形成了一些强疏水性接触,尤其是 Fi19-Fj19和Fi19-Fj20,其中 i 和 j 指的是寡聚体中的两个不同肽链。

图 2  6 个 Aβ16-22 肽聚集后的 100 ns MD 模拟分析

(A)显示了系统随时间的低聚状态;(B)具有概率的残基间接触图。

作者: 蒲中机

中机是大连理工大学在读博士,擅长生物体系分子模拟领域,精通软件Amber,Gromacs,autodock,ledock,modeller ...;目前运营微信公众号生物研究生,搜索公众号shengwuyanjiusheng,可关注



发表回复

您的电子邮箱地址不会被公开。

联系我们

联系我们

(44)07934433023

在线咨询: QQ交谈

邮箱: info@bioengx.org

关注微信
微信扫一扫关注我们

微信扫一扫关注我们

关注微博
返回顶部
Designed by

best down free | web phu nu so | toc dep 2017