本实验方案采用一种生物分子模拟软件，描述了用于模拟野生型半胱天冬酶（caspase）及其突变体的分子动力学（MD）方法。分子动力学方法可用于评估半胱天冬酶结构的动态演化过程，以及突变或翻译后修饰可能产生的影响。
本文内容
摘要
摘要
细胞凋亡是一种程序性细胞死亡，可清除受损细胞，并调控多细胞生物的发育过程和组织稳态。半胱天冬酶（caspase）是一类半胱氨酸蛋白酶家族，在细胞凋亡的启动和执行过程中发挥关键作用。半胱天冬酶的成熟及其活性通过翻译后修饰以高度动态的方式进行精细调控。为了评估翻译后修饰的影响，通常将潜在修饰位点突变为不受任何修饰影响的氨基酸残基。例如，将丝氨酸残基替换为丙氨酸或天冬氨酸。然而，此类替换可能改变半胱天冬酶活性位点的构象，从而干扰其催化活性和细胞功能。此外，位于关键位置的其他氨基酸残基的突变也可能破坏半胱天冬酶的结构与功能，导致细胞凋亡过程紊乱。为避免使用突变残基所带来的困难，可直接采用分子建模方法来估算氨基酸替换对半胱天冬酶结构的潜在影响。本实验方案利用生物分子模拟软件包（Amber）和超级计算机资源，对野生型半胱天冬酶及其突变体进行建模，以检测突变对蛋白质结构和功能的影响。
引言
细胞凋亡是调控多细胞生物形态发生和组织稳态研究最为广泛的一种细胞过程。细胞凋亡可由多种外部或内部刺激引发，例如死亡受体的激活、细胞周期信号的紊乱、DNA损伤、内质网（ER）应激以及多种细菌和病毒感染1。半胱天冬酶（caspases）是细胞凋亡中的关键执行分子，通常根据其结构域特征及在caspase级联反应中的位置分为两类：起始型（initiator caspases，包括caspase-2、caspase-8、caspase-9和caspase-10）和效应型（effector caspases，包括caspase-3、caspase-6和caspase-7）2,3。当接收到细胞死亡信号时，起始型caspase会与衔接分子相互作用，促进其近端诱导的二聚化和自剪切，从而形成有活性的酶。效应型caspase则通过起始型caspase的切割而被激活，并通过裂解多种细胞底物来执行下游的凋亡程序4。
起始型和效应型caspase的成熟与功能受到多种不同细胞内机制的调控，其中翻译后修饰在细胞死亡调控中发挥着不可或缺的作用5。添加修饰基团（如磷酸化、亚硝基化、甲基化或乙酰化）或蛋白质（如泛素化或SUMO化）可改变caspase的酶活性，或影响其蛋白构象与稳定性，从而调控凋亡过程。定点突变技术被广泛用于研究潜在的翻译后修饰位点并阐明其功能。通常将假定的修饰位点突变为不能被进一步修饰的其他氨基酸，例如将可能被磷酸化的丝氨酸和苏氨酸突变为丙氨酸，或将赖氨酸的泛素化位点替换为精氨酸。另一种策略是引入特定氨基酸以模拟翻译后修饰状态（例如，谷氨酸和天冬氨酸已被用于模拟磷酸化的丝氨酸或苏氨酸）6。然而，某些位于活性位点附近或关键位置的氨基酸替换可能改变caspase的结构，干扰其催化活性，并抑制凋亡性细胞死亡7。在caspase基因的肿瘤相关错义突变中也可观察到类似效应。例如，caspase-6的肿瘤相关突变R259H导致底物结合口袋中环状结构的构象改变，降低了底物的有效催化转换率8。在头颈部鳞状细胞癌中发现的caspase-8的G325A氨基酸替换可能损害caspase-8的活性，进而调控核因子-kB（NF-kB）信号通路，促进肿瘤发生9。
为了评估氨基酸替换对半胱天冬酶（caspase）结构和功能的潜在影响，可采用分子建模方法。本文描述了使用生物分子模拟软件包（Amber）对野生型半胱天冬酶及其突变体进行建模的分子动力学（MD）方法。分子动力学方法能够展示引入突变后蛋白质结构的动态演化过程。Amber软件包最初由Peter Kollman团队开发，现已成为生物分子模拟中最常用的软件工具之一10,11,12,13。该软件分为两个部分：（1）AmberTools，是一组常用于系统准备（如原子类型分配、添加氢原子和显式水分子等）和轨迹分析的程序；（2）Amber，其核心是pmemd模拟程序。AmberTools为免费软件包（同时也是安装Amber的前提条件），而Amber则需通过独立的许可协议并支付费用获取。在超级计算机上或利用图形处理器（GPU）进行并行模拟，可显著提升蛋白质结构动力学科学研究的计算性能14。目前最新的可用软件版本为AmberTools21和Amber20，但所描述的实验方案同样适用于早期版本。
方案
1. 系统准备
注意：天然蛋白和突变体蛋白的分子模型是基于从蛋白质数据库（Protein Data Bank）获得的合适晶体结构构建的15,16。
- 为获取选定的 PDB 结构，使用 Download Files 下拉菜单并点击 PDB Format。从 PDB 文件中移除注释和连接性数据，并在不同蛋白质链之间插入 TER 卡片。可选地，通过将组氨酸残基名称 HIS 替换为 HIE、HID 或 HIP 来设定其质子化状态（补充文件 1）。
注意：保留晶体结构中解析出的水分子（如果 PDB 文件中存在）是合理的。 - 为准备起始模型，启动 tleap 程序（AmberTools 软件包），然后输入用于操作 tleap 对象的命令。
- 启动程序：tleap。加载 ff14SB 力场以用分子力学方法描述蛋白质：> source leaprc.protein.ff14SB。
- 加载水分子及原子离子（Na+、Cl-）的力场：> source leaprc.water.tip3p。
- 加载 PDB 文件 并构建氢原子坐标，创建一个名为 mol 的对象：> mol = loadpdb protein.pdb。
- 检查可能引起问题的内部不一致性：> check mol。
注意：如果存在二硫键，必须手动指定。将共价连接的半胱氨酸残基名称 CYS 替换为 CYX，并在 tleap 中输入以下命令以在 SG 原子间建立化学键：bond mol.X.SG mol.Y.SG（X 和 Y 为半胱氨酸残基编号）。 - 在蛋白质周围创建溶剂盒（TIP3P 水模型，12 Å 距离）：> solvatebox mol TIP3PBOX 12.0。
- 检查总电荷：> charge mol，并添加抗衡离子（Na+ 或 Cl-）以中和体系：
> addions mol Na+ 0。 - 生成拓扑 prmtop 文件 和坐标 inpcrd 文件（用于运行模拟的输入文件）：> saveamberparm mol protein.prmtop protein.inpcrd。退出 tleap 程序：> quit。
注意：可使用 ambpdb 工具轻松将 Amber 格式的坐标文件转换为 PDB 格式（参见《用户手册》，见 材料表）。
2. 能量最小化
注意：能量最小化是必要的，可消除初始体系中原子间的不良接触和重叠，这些因素在进行分子动力学模拟时会导致体系不稳定。
- 执行能量最小化的第一阶段（最陡下降法 2,500 步 + 共轭梯度法 2,500 步），在对重原子施加位置限制以固定蛋白质坐标的同时，优化所添加的氢原子和水分子的位置。
- 按如下方式运行 pmemd 程序：
pmemd -O -i min1.in -p protein.prmtop -c protein.inpcrd -o protein_min1.out
-r protein_min1.rst -ref protein.inpcrd。 - 遵循所需的参数：-i FILE 控制数据；-p FILE 分子拓扑、力场参数、原子名称；-c FILE 初始坐标；-o FILE 用户可读的日志输出；-r FILE 最终坐标；-ref FILE 用于位置限制的参考坐标。
- 在输入文件 min1.in 中指定选项：
&cntrl
imin=1, maxcyc=5000, ncyc=2500,
cut=10.0, ntb=1,
ntc=1, ntf=1,
ntpr=10,
ntr=1,
restraintmask=':1-517 & !@H=',
restraint_wt=2.0
/
注意：输入选项说明：imin=1 执行最小化；maxcyc=5000 最小化循环的最大次数；ncyc=2500 在 ncyc 次循环后，最小化方法从最陡下降法切换为共轭梯度法；cut=10.0 非键截断距离（Å）；ntb=1 施加周期性边界条件，恒定体积；ntc=1 对键长不施加约束（以获得更好的能量收敛性）；ntf=1 计算所有相互作用；ntpr=10 每 ntpr 步将能量信息输出至 out 文件；ntr=1 使用谐波势能限制指定原子的标志；restraintmask=':1-517 & !@H=' 指定受限制的原子；restraint_wt=2.0 位置限制的权重（kcal/mol∙Å2）。
- 按如下方式运行 pmemd 程序：
- 执行能量最小化的第二阶段（最陡下降法 5,000 步 + 共轭梯度法 5,000 步），在无限制条件下优化整个体系。
- 使用 min2.in 和 protein_min1.rst 作为输入文件：pmemd -O -i min2.in -p protein.prmtop
-c protein_min1.rst -o protein_min2.out -r protein_min2.rst。 - 在输入文件 min2.in 中指定选项：
&cntrl
imin=1, maxcyc=10000, ncyc=5000,
cut=10.0, ntb=1,
ntc=1, ntf=1,
ntpr=10
/
- 使用 min2.in 和 protein_min1.rst 作为输入文件：pmemd -O -i min2.in -p protein.prmtop
3. 加热
注意：此阶段旨在将体系从 0 K 加热至 300 K。由于基于 PDB 文件的初始模型不包含原子速度信息，因此需要为原子分配初始速度。
- 对蛋白质原子施加位置限制进行加热过程（50 ps，恒定体积）。
- 使用 heat.in 和 protein_min2.rst 作为输入文件：pmemd -O -i heat.in -p protein.prmtop
-c protein_min2.rst -o protein_heat.out
-r protein_heat.rst -x protein_heat.mdcrd -ref protein_min2.rst - 遵循所需参数：-x FILE 表示在分子动力学轨迹中保存的坐标集合。
- 在输入文件 heat.in 中指定选项：
&cntrl
imin=0, irest=0, ntx=1,
nstlim=25000, dt=0.002,
ntc=2, ntf=2,
cut=10.0, ntb=1,
ntpr=500, ntwx=500,
ntt=3, gamma_ln=2.0,
tempi=0.0, temp0=300.0,
ntr=1, restraintmask=':1-517',
restraint_wt=1.0,
nmropt=1
/
&wt TYPE='TEMP0',
istep1=0, istep2=25000,
value1=0.1, value2=300.0 /
&wt TYPE='END' /
注：输入参数说明：imin=0 执行分子动力学模拟；irest=0 开始新的模拟；ntx=1 不从 rst 文件中读取初始速度；nstlim=25000 分子动力学步数；dt=0.002 时间步长（ps）；ntc=2 使用 SHAKE 算法约束含氢键的键长；ntf=2 不计算被约束键的力；ntwx=500 每 ntwx 步将坐标写入 mdcrd 文件；ntt=3 使用 Langevin 恒温器；gamma_ln=2.0 Langevin 动力学的碰撞频率（ps−1）；tempi=0.0 初始温度；temp0=300.0 参考温度；nmropt=1 读取 &wt 名单中指定的参数；TYPE='TEMP0' 目标温度逐步变化。
- 使用 heat.in 和 protein_min2.rst 作为输入文件：pmemd -O -i heat.in -p protein.prmtop
4. 平衡
注意：此步骤对于调节水的密度并获得蛋白质的平衡状态是必要的。
- 在300 K下进行平衡，不施加任何限制（500 ps，恒压）。
- 使用 平衡液 和 protein_heat.rst 作为输入： pmemd -O -i 平衡.in -p 蛋白质.prmtop
-c protein_heat.rst -o protein_equil.out -r protein_equil.rst -x protein_equil.mdcrd. - 在输入文件中指定选项 平衡中:
&cntrl
imin=0，irest=1，ntx=5，
nstlim=250000,
dt=0.002,
ntc=2，ntf=2，
cut=10.0，ntb=2，ntp=1，taup=2.0
ntpr=1000，ntwx=1000，ntwr=50000，
ntt=3，gamma_ln=2.0，
temp0=300.0
/
注意：输入选项： irest=1 重新启动模拟； ntx=5 坐标和速度从先前生成的文件中读取 rst 文件； ntb=2 周期性边界 施加恒定压力； ntp=1 各向同性压力缩放； τau_p = 2.0 压力弛豫时间（ps）； ntwr=50000 每个 ntwr 步骤，这些 rst 文件已写入，确保从崩溃中恢复。
- 使用 平衡液 和 protein_heat.rst 作为输入： pmemd -O -i 平衡.in -p 蛋白质.prmtop
5. 产量动态
- 达到平衡后，进行恒压条件下的生产阶段分子动力学模拟（10 ns 或更长时间），并生成轨迹文件，用于后续分析蛋白质结构。
- 使用 prod.in 和 protein_equil.rst 作为输入文件：pmemd -O -i prod.in -p protein.prmtop
-c protein_equil.rst -o protein_prod.out -r protein_prod.rst -x protein_prod.mdcrd。 - 在输入文件 prod.in 中指定选项：
&cntrl
imin=0, irest=1, ntx=5,
nstlim=5000000,
dt=0.002,
ntc=2, ntf=2,
cut=10.0, ntb=2, ntp=1, taup=2.0,
ntpr=1000, ntwx=1000, ntwr=50000,
ntt=3, gamma_ln=2.0,
temp0=300.0, ig=-1
/
注意：可在计算机集群和超级计算机上使用 pmemd 的并行版本（pmemd.MPI）或 GPU 加速版本（pmemd.cuda）。长时间的分子动力学模拟可拆分为多个阶段依次执行。强烈建议每次重启模拟时设置 ig=-1（随机种子选项）。可使用 ptraj 或 cpptraj 程序对坐标轨迹进行分析和处理。更多详细信息，请参阅软件用户手册。
- 使用 prod.in 和 protein_equil.rst 作为输入文件：pmemd -O -i prod.in -p protein.prmtop
结果
本方案可直接应用于半胱天冬酶（caspase）的翻译后修饰或致病性突变的研究。本节展示了分子动力学（MD）建模的工作流程（图1），该方法已成功用于caspase-2的研究7。通过体外（in vitro）定点突变技术将潜在的磷酸化位点（Ser/Thr突变为Ala）并结合生物化学方法，研究发现Ser384Ala突变可阻止caspase-2的加工过程，并抑制其酶活性以及凋亡性细胞死亡的诱导（补充图1）。然而，质谱分析和Phos-Tag技术均未证实Ser384位点发生磷酸化7。为了阐明Ser384如何调控caspase-2的活性，研究人员对该蛋白的三维结构进行了分子动力学（MD）建模分析（图1）。
使用1pyo晶体结构构建了野生型和突变型caspase-2的分子模型（可用的caspase-2结构列于补充表1中）。Ser384Ala突变体通过去除Ser384残基中的Oγ原子而构建（在PDB中，Ser和Ala残基的区别在于是否含有Oγ原子）。晶体结构中通常不包含氢原子坐标，因此向蛋白质结构中添加了氢原子，随后用厚度为12 Å的TIP3P水层进行溶剂化，以进行后续的显式溶剂模拟。添加钠离子以中和体系电荷。所得的caspase-2起始模型根据分子动力学（MD）建模流程，使用min1.in、min2.in、heat.in、equil.in和prod.in控制数据文件，依次进行能量最小化、平衡和随后的10 ns MD模拟。
在caspase-2的晶体结构中，Ser384与Arg219和Arg378共同参与构成活性位点腔体的表面。Arg219和Arg378与底物的羧基形成氢键，而Ser384不直接与底物发生相互作用。分子动力学（MD）模拟证实，Ser384Ala取代未影响催化残基Cys320（亲核体）和His277（广义碱），但导致Arg378发生显著的构象变化。其胍基转向体相溶剂，使得Nε原子无法与底物的羧基形成关键的氢键（图1）。在天然酶中，Ser384的Oγ原子（部分负电荷）与Arg378的胍基（正电荷）之间存在静电相互作用。然而，丝氨酸的取代显然破坏了这一相互作用。因此，研究表明Ser384Ala取代影响了caspase-2活性位点中精氨酸残基对底物的识别，从而削弱了酶活性以及诱导凋亡性细胞死亡的能力。该发现的机制在进化上似乎保守，并在其他caspase家族成员中普遍存在。因此，分子动力学模拟的应用不仅验证了生化实验结果，还进一步揭示了caspase活性中心分子结构的新见解。
图1：用于研究caspase结构的分子动力学建模流程。 依据方案部分的说明，对野生型caspase-2及其Ser384Ala突变体进行了研究。结果显示，Ser384Ala替换在活性位点残基Arg378处诱导了重要的构象变化。 请点击此处查看该图的放大版本。
补充图1：caspase-2在细胞死亡中的功能。 在外源性和内源性刺激作用下，野生型caspase-2可通过自蛋白水解机制被激活。活化的caspase-2切割Bid蛋白，进而促进线粒体外膜通透化，诱导细胞凋亡。Ser384Ala突变可阻止caspase-2的活化及其介导的细胞死亡。 请点击此处下载该文件。
补充表 1：蛋白质数据库中可获取的 Caspase-2 结构。 结构按发布日期排序。请点击此处下载该表格。
补充文件1：编辑PDB文件的不同步骤。 请点击此处下载该文件。
讨论
所述的分子动力学方法可用于对胱天蛋白酶（caspase）的野生型和突变型进行建模，借助生物分子模拟软件包实现。本文讨论了该方法中若干重要的技术问题。首先，需从蛋白质数据库（Protein Data Bank）中选择一个具有代表性的胱天蛋白酶晶体结构。值得注意的是，胱天蛋白酶的单体和二聚体形式均可接受。建议选择分辨率较高且缺失残基数量最少的结构。某些残基的质子化状态可在输入的PDB文件中手动设定。例如，在图1中，催化残基His277的Nε2原子上连接了一个氢原子（HIE形式）。其次，必须在蛋白质周围构建溶剂盒子，以进行后续的显式溶剂模拟。需要进行能量最小化，以优化所添加的氢原子和水分子的坐标。第三，在开始生产性动力学模拟之前，必须完成加热和平衡过程。在加热阶段，应确保蛋白质表面的空腔和结合口袋被水分子充分填充（如有必要，可增加分子动力学步数）。在平衡阶段，需通过分析蛋白质主链原子相对于初始位置的均方根偏差来确认其达到平衡构象12。第四，在生产性动力学模拟过程中，必须监测蛋白质结构的构象变化。为此，应通过拟合主链原子将轨迹帧与起始结构叠加，然后使用分子可视化工具（如VMD、PyMOL等）分析构象变化17。如有需要，可增加分子动力学步数，或将其拆分为多个连续执行的模拟段。第五，建议使用GPU加速版本的分子动力学模拟程序（pmemd.cuda），其性能显著优于传统的CPU实现方式14,18。
所述方法的一个可能局限性在于蛋白质数据库（Protein Data Bank）中某些caspase结构的可获得性。然而，在PDB中已有超过900条条目属于该类别 "caspase" 请求表明，caspase 的结构已被深入研究。当晶体结构的分辨率较低（即细节层次粗糙）或缺少对 caspase 功能至关重要的某些片段时，可能会出现一些困难。
总之，分子动力学（MD）建模方法已被证明在预测基因突变、翻译后修饰或分子间相互作用后蛋白质结构变化方面具有有效性。在细胞死亡领域应用MD模拟，为研究翻译后修饰调控半胱天冬酶（caspase）的分子机制提供了重要机遇。在此背景下，长时间的MD模拟可用于评估功能区域的动态行为以及氨基酸替换的潜在影响，从而阐明与突变或修饰的caspase-2及其他caspase相关的分子机制。例如，在胃肠道癌症以及中枢和周围神经系统肿瘤中已检测到起始型caspase基因（caspase-2/caspase-8/caspase-9/caspase-10）的错义突变19。综上所述，对caspase进行MD建模是一种强大的in silico方法，有助于理解调控程序性细胞死亡及其相关疾病发生的机制，并推动新型有效疗法的开发。
披露
作者声明无利益冲突。
致谢
本研究得到了俄罗斯科学基金会（项目编号：17-75-20102，用于实验方案开发）的资助。代表性结果部分所述实验（磷酸化分析）得到了斯德哥尔摩癌症协会（181301）和瑞典癌症协会（190345）的支持。
材料
| 姓名 | 公司 | 目录编号 | 评论 |
|---|---|---|---|
| Amber20 | 加利福尼亚大学旧金山分校 | 用于分子动力学模拟的软件 http://ambermd.org | |
| AmberTools21 | 加利福尼亚大学旧金山分校 | 用于分子建模与分析的软件 http://ambermd.org |
参考文献
- Olsson, M., Zhivotovsky, B. Caspases and cancer. Cell Death and Differentiation. 18 (9), 1441-1449 (2011).
- Lavrik, I. N., Golks, A., Krammer, P. H. Caspase: Pharmacological manipulation of cell death. Journal of Clinical Investigation. 115 (10), 2665-2672 (2005).
- Degterev, A., Boyce, M., Yuan, J. A decade of caspases. Oncogene. 22 (53), 8543-8567 (2003).
- Pop, C., Salvesen, G. S. Human caspases: Activation, specificity, and regulation. Journal of Biological Chemistry. 284 (33), 21777-21781 (2009).
- Zamaraev, A. V., Kopeina, G. S., Prokhorova, E. A., Zhivotovsky, B., Lavrik, I. N. Post-translational modification of caspases: The other side of apoptosis regulation. Trends in Cell Biology. 27 (5), 322-339 (2017).
- Pearlman, S. M., Serber, Z., Ferrell, J. E. A mechanism for the evolution of phosphorylation sites. Cell. 147 (4), 934-946 (2011).
- Zamaraev, A. V., et al. Requirement for Serine-384 in Caspase-2 processing and activity. Cell Death and Disease. 11 (10), 825(2020).
- Dagbay, K. B., Hill, M. E., Barrett, E., Hardy, J. A. Tumor-associated mutations in caspase-6 negatively impact catalytic efficiency. Biochemistry. 56 (34), 4568-4577 (2017).
- Ando, M., et al. Cancer-associated missense mutations of caspase-8 activate nuclear factor-κB signaling. Cancer Science. 104 (8), 1002-1008 (2013).
- Case, D. A., et al. AMBER 2020. , University of California, San Francisco. (2020).
- Salomon-Ferrer, R., Case, D. A., Walker, R. C. An overview of the Amber biomolecular simulation package. WIREs Computational Molecular Science. 3 (2), 198-210 (2013).
- Roe, D. R., Cheatham, T. E. PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data. Journal of Chemical Theory and Computation. 9 (7), 3084-3095 (2013).
- Maier, J. A., et al. ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB. Journal of Chemical Theory and Computation. 11 (8), 3696-3713 (2015).
- Salomon-Ferrer, R., Götz, A. W., Poole, D., Le Grand, S., Walker, R. C. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh ewald. Journal of Chemical Theory and Computation. 9 (9), 3878-3888 (2013).
- Berman, H. M., et al. The Protein Data Bank. Nucleic Acids Research. 28 (1), 235-242 (2000).
- Burley, S. K., et al. RCSB Protein Data Bank: Powerful new tools for exploring 3D structures of biological macromolecules for basic and applied research and education in fundamental biology, biomedicine, biotechnology, bioengineering and energy sciences. Nucleic Acids Research. 49, 437-451 (2021).
- Martinez, X., et al. Molecular graphics: Bridging structural biologists and computer scientists. Structure. 27 (11), 1617-1623 (2019).
- Lee, T. S., et al. GPU-accelerated molecular dynamics and free energy methods in Amber18: Performance enhancements and new features. Journal of Chemical Information and Modeling. 58 (10), 2043-2050 (2018).
- Jäger, R., Zwacka, R. M. The enigmatic roles of caspases in tumor development. Cancers. 2 (4), 1952-1979 (2010).
重印与许可
