使用MolAICal基于NAMD模拟结果计算小分子和蛋白MM/PBSA的教程

MolAICal

之前发帖是MM/GBSA: http://bbs.keinsci.com/thread-18890-1-1.html , 这一版是MM/PBSA的计算教程。

目前，还可以基于单个pdb文件计算MM/PBSA或MM/GBSA；由于发帖字数限制-->本教程省略了部分参数，更详细的教程见：https://molaical.gitlab.io/cntutorial.html  或  https://molaical.github.io/tutorial.html

1. 简介

在本教程中介绍了基于NAMD的分子动力学模拟结果，使用MolAICal计算小分子和Mpro蛋白受体MM/PBSA的方法。本教程只是一个简单演示。为了节省运行及存储空间，本教程仅选择了Mpro复合物分子动力学模拟的25帧用于计算。本教程不仅可以用于计算蛋白质-配体的MM/PBSA，还可以用于计算基于MD模拟的蛋白质-多肽、蛋白质-蛋白质、DNA-配体、DNA-蛋白质、蛋白质-DNA、蛋白质-RNA和其它任意复合物的MM/PBSA；只需要用指定的对象替换本教程中的蛋白质和配体即可。例如，基于本教程中相同的运行命令参数，本教程中的蛋白质被替换为DNA，配体被替换为多肽。

2.工具

2.1. 所需软件下载地址

1) MolAICal: https://molaical.github.io  或  https://molaical.gitee.io

2) NAMD: https://www.ks.uiuc.edu/Research/namd/

注意：本教程可以使用NAMD2.x、3.x或更高版本。例如，如果使用NAMD 3.x版本，则使用命令“namd3”替代本教程中的命令“namd2”。对于更高版本的NAMD，用户可以使用与前面示例类似的替换方式。

3) Carma: https://github.com/glykos/carma or https://utopia.duth.gr/~glykos/Carma.html

4) APBS (https://apbs.readthedocs.io) or DelPhi (http://compbio.clemson.edu/lab/delphisw)

2.2. 操作示例文件

所有用到的操作教程文件均可在下面的网站下载：

https://gitee.com/molaical/tutorials/tree/master/020-mmpbsa

3. 操作流程

第一部分：问题的解决

问题: 有些用户使用MM/PBSA进行计算得到正值，如下:

https://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2020-2021/1295.html

如果用户的研究对象不存在多个分子进入不同镜像的问题，可以跳过这个部分，直接进入Part II.

可能的解决方案：必须检查一下配体是否一直跟受体在同一个周期性盒子中，如果不是，请使用VMD PBC将配体和受体裹进同一个周期性盒子中。可以参考：

https://www.ks.uiuc.edu/Research/vmd/script_library/scripts/pbcwrap/

https://www.ks.uiuc.edu/Research/vmd/plugins/pbctools/

这里我提供2个命令在VMD Tk Console中运行:

1. %> pbc wrap -centersel protein -center com -compound chain -all

复制代码

注：在VMD中检查复合物结构，若上述命令能将复合物整体包裹在一起（即：去周期化条件），则可跳过以下命令（pbc unwrap）。

1. %> pbc unwrap -sel "not (water or ions)" -all

复制代码

运行上述命令可以将配体和受体包进同一个周期性盒子。使用下面的命令保存轨迹到DCD文件中：

1. %> set all [atomselect top all]
   
2. %> animate write dcd pro.dcd sel $all beg 0 end 24 waitfor all

复制代码

这里所有原子都被选中，用户可以根据需求进行选择。为了简便，本教程只提供25帧：从0到24帧。"pro.dcd" 是轨迹名字。解决上述问题之后，可以开始下面的教程。

选项（无图形界面的去周期化条件操作）

若用户已完成去周期化条件操作，或已经在Linux版本的MolAICal容器中配置好VMD和NAMD，或无需对DCD文件进行去周期化条件操作，可直接跳过此选项步骤。

用户可在 Windows 或 Linux 的图形桌面环境中执行并测试上述命令。但在无图形界面的环境中操作时，需通过命令行输入。为支持外部程序调用，MolAICal 提供了持久化配置接口。以下命令用于配置 VMD（Visual Molecular Dynamics）的名称与路径，仅需单次配置即可供后续重复使用，无需重复设置。

为了解决Linux中的库依赖问题，Linux版本的MolAICal (Windows 版本的MolAICal没有采用容器)采用了基于容器的方法。如果需要调用外部程序，建议在MolAICal容器内安装这些程序。如果不需要外部程序，可以忽略此步骤。本教程需要外部程序NAMD和VMD，具体步骤如下：

1. 首先，将文件复制到 MolAICal 容器中

1. # 进入容器文件系统（将进入 "/root" 目录）
   
2. #> molaical.exe -eset shell in
   
3. 
   
4. # 将 VMD 和 NAMD 安装包从本地机器复制到容器中，'cp' 命令的第一部分（源路径）位于本地主机，第二部分（目标路径）在容器内；VMD 和 NAMD 软件包可通过以下命令移入容器。
   
5. #> cp /home/user/<本地文件> /root/soft
   
6. 
   
7. # 退出容器文件系统
   
8. #> exit

复制代码

2.      其次，进入 MolAICal 容器的虚拟环境

1. #> molaical.exe -eset sys run molaical

复制代码

注：molaical 是容器名称。

3. 在 MolAICal 容器虚拟环境中安装软件的方式与在本地主机上安装相同。以下以安装 VMD 和 NAMD 为例：

1) 安装NAMD：

解压 NAMD 文件（假设解压后的文件夹名为 namdcpu），然后使用以下命令将其路径告知 MolAICal：

1. #> molaical.exe -call set -n NAMD -p "/root/soft/namdcpu/namd3"

复制代码

注：请将上述 VMD 和 NAMD 的路径替换为您系统中的实际路径。-n 后面的 "VMD" 和 "NAMD"（大小写不敏感）是固定的标识符。

2) 安装VMD：

按以下步骤操作：

• 解压 VMD 文件：
  

1. #> tar -xzvf vmd-xxx.tar.gz

复制代码

注：请将上述路径替换为您系统中的实际路径。

• 修改 VMD 解压目录中名为configure的文件中的安装路径：
  

1. # 默认值：
   
2. $install_bin_dir="/usr/local/bin";
   
3. $install_library_dir="/usr/local/lib/$install_name";
   
4. 
   
5. # 修改为：
   
6. $install_bin_dir="/root/soft/vmd193/bin";
   
7. $install_library_dir="/root/soft/vmd193/lib/$install_name";

复制代码

• 安装 VMD：
  

1. #> cd vmd-xxx
   
2. #> ./configure LINUXAMD64
   
3. #> cd src
   
4. #> make install

复制代码

注：请运行 ./configure 并根据所用计算机选择正确的类型，此处为 "LINUXAMD64"。

• 然后使用以下命令将 VMD 路径告知 MolAICal：
  

1. #> molaical.exe -call set -n VMD -p "/root/soft/vmd193/bin/vmd"

复制代码

至此，NAMD 和 VMD 在 MolAICal 容器内的安装与配置已完成。为防止MolAICal出现问题时丢失已安装的程序（例如VMD和NAMD），请参阅附录1中的第3步。

4. 记得使用exit命令退出MolAICal虚拟环境，返回本地计算机进行计算（主要是为了省去拷贝文件的步骤，在容器内运行也行，但需要从本地计算机向MolAICal容器中拷贝数据。）

1. #> exit

复制代码

然后，运行外部程序的参数：

1. #> molaical.exe -call run -n vmd -i -dispdev text -e  pbcwrap.vmd

复制代码

注意：

1) 在 MolAICal 的环境变量配置中，'-n' 后的参数名称（如 '-n VMD'）必须与对应的运行时参数（如 '-n vmd'）拼写一致，但大小写可忽略。如上述两个命令中'-n' 后的参数名称需保持一致。

2) '-i' 后接 VMD 的运行时参数，且 '-i' 必须置于其他参数之后。此处，用户只需将 pbcwrap.vmd 替换为自己的脚本文件，当然也可直接选用本教程提供的脚本 pbcwrap.vmd。

下方以 pbcwrap.vmd 为例（用户可根据实际用途或分析结果修改此脚本）：

1. package require pbctools
   
2. mol new mpro.psf
   
3. mol addfile mpro.dcd waitfor all
   
4. pbc wrap -centersel protein -center com -compound chain -all
   
5. set all [atomselect top all]
   
6. animate write dcd mpro.dcd sel $all beg 0 end 24 waitfor all
   
7. quit

复制代码

第二部分：基于NAMD轨迹进行MM/PBSA的计算

假设用户已在MolAICal中配置了VMD和NAMD路径。上述说明适用于Linux版MolAICal容器的外部程序安装。对于Windows版MolAICal，用户需将以下命令中的路径替换为实际系统路径后直接运行：

1. #> molaical.exe  -call set -n VMD -p "C:\Program Files (x86)\University of  Illinois\VMD\vmd.exe"  
   
2. #> molaical.exe  -call set -n namd -p "d:\namd2\namd2.exe"  
   
3. #> molaical.exe  -call set -n delphi -p "E:/delphi/delphicpp_v8.4.3_windows_x64.exe"

复制代码

注意：请将上述 VMD 和 NAMD 的路径替换为用户系统中的实际路径。当运行MMPB/GB/SA功能并通过MolAICal设置环境变量时（例如：#> molaical.exe -call set -n VMD -p"path"），MMPB/GB/SA中与"-n"参数对应的变量名是固定的，包括：'vmd'、'namd'和'delphi'（不区分大小写）。

• 其中DelPhi是二进制文件，在linux系统上需要赋予可执行权限： chmod +xdelphicpp_v8.4.3_windows_x64.exe ，所以，可以参考附录1中的NAMD安装步骤。
• 请注意，Windows版与Linux版的MolAICal配置存在差异：Linux版本采用基于udocker容器的技术方案，需在容器内部完成设置，而Windows版本则无需此步骤。
  

3.1 MM/PBSA的计算

假设已经有跑平衡的轨迹名为"mpro.dcd"，同时，在"mpro.dcd"轨迹中，保证受体和配体在同一个周期性盒子中。用户可以替换成自己的轨迹。转到以下目录：

1. #> cd 020-mmpbsa

复制代码

然后运行以下命令，可以快速计算MM/PBSA:

1. #> molaical.exe -mmpbgbsa -fmmpbsa_apbs.vmd

复制代码

• 可通过修改文件'mmpbsa_apbs.vmd'中的参数，或通过输入'-i'参数进行修改。如需更详细信息，请查阅MolAICal手册。
• 选项：此处输入文件'mmpbsa_apbs.vmd'中的'mpro.dcd'和'mpro.psf'等文件包含蛋白质、配体、离子、水分子等成分，可能会占用大量计算机内存。为节省计算内存，用户可使用以下命令仅提取受体与配体复合物的轨迹：
  

1. #> molaical.exe -call run -n vmd -c stripDCD -i -dispdev  text -psf "mpro.psf" -args protein,or,resname,LIG  "mpro.dcd" "complex" mpro.psf mpro.pdb

复制代码

注："complex"是输出文件的前缀。系统将参考mpro.psf和mpro.pdb文件生成以"complex"为前缀的文件，包括：'complex.psf'、'complex.pdb'和'complex.dcd'。这些文件可替代输入文件'mmpbsa_apbs.vmd'中的'mpro.psf'、'mpro.pdb'和'mpro.dcd'，从而节省加载内存。

运行上述命令后，MM/PBSA的输出结果和结合自由能△G如下:

------------------------------------------------------------------------------------------------------------------

   Pol:                           25               -6.6470                  2.6141         13.0705

   Npol:                         25               -59.8536                1.1223         5.6116

   G binding:                25               -66.5006        +/-  2.7810          13.9052

-----------------------------------------------------------------------------------------------------------------

注：所有能量值单位均为千卡/摩尔（kcal/mol）。在不同操作系统或NAMD版本中计算结果可能存在差异，但若保留两位小数，结果将保持一致或高度吻合。

输入文件与结果分析

打开 ‌'mmpbsa_apbs.vmd'‌ 文件，请注意这两个参数：‌"-poisson_boltzmann"和‌"-pb_executable":‌

• 若用户未设置 ‌"-pb_executable"‌ 参数，MolAICal 将自动调用内置的 APBS 计算泊松-玻尔兹曼（PB）能量。
• 若 ‌"-poisson_boltzmann"‌ 设为 ‌2‌，程序将调用 APBS 计算 PB 能量；设为 ‌1‌ 则调用 DelPhi 计算 PB 能量；设为 ‌0‌ 则执行 MM/GBSA 计算。
  

若用户将 ‌"-poisson_boltzmann"‌ 设为 ‌1‌，并将 ‌"-pb_executable"‌ 指向 DelPhi 的路径（本教程中文件名为 ‌'mmpbsa_delphi.vmd'‌），再次运行以下命令时，程序将调用 DelPhi 计算 PB 能量：

1. #> molaical.exe -mmpbgbsa -f mmpbsa_delphi.vmd

复制代码

• 其中，文件'mmpbsa_apbs_create_file.vmd'‌或‌'mmpbsa_delphi_create_file.vmd'‌比 ‌'mmpbsa_apbs.vmd'‌ 或 ‌'mmpbsa_delphi.vmd'‌ 多出以下参数：
• "-namd_config_file      input.namd \".
• "-pb_gb_file             input.apbs \"   或"-pb_gb_file            input.delphi \"
• "-complex_selection"‌、‌"-receptor_selection"‌ 和 ‌"-ligand_selection"‌ 用于选择复合物、受体和配体的分子部分，功能类似 VMD Tk 控制台中的 ‌"atomselect"‌ 命令。
• "-sa_input_apbs"‌ 用于自定义 SASA 计算的 APBS 输入文件。用户可参考材料文件中的模板 ‌"input_sa.apbs"‌，只需在输入文件（如 ‌'mmpbsa_apbs.vmd'‌ 或 ‌'mmpbsa_delphi.vmd'‌）中添加一行："-sa_input_apbs      input_sa.apbs \"
  

通过此机制，用户可以在 ‌"input.namd"‌ 文件中自定义 NAMD 的输入参数，在 ‌"input.apbs"‌ 文件中自定义 APBS 的输入参数，或在 ‌"input.delphi"‌ 文件中自定义 DelPhi 的输入参数。若运行以下命令：

1. #> molaical.exe -mmpbgbsa -f mmpbsa_apbs_create_file.vmd
   
2. 或  
   
3. #> molaical.exe -mmpbgbsa -f mmpbsa_delphi_create_file.vmd

复制代码

且自定义文件（如："input.namd"或"input.apbs"）中的参数与 MolAICal 内置参数一致，则将得到相同结果。

• 如上表结果所示，极性贡献为 ‌Pol = Elec + PB‌，非极性贡献为 ‌Npol = Vdw + SA‌。
  

‌

第三部分：残基能量分解

3.2. 残基能量分解

可以参考上面的方法进行残基能量分解，比如：根据报道，SARS-CoV-2 Mpro的残基M165和配体N3有相互作用，本教程就以残基M165为例，计算残基M165的自由能贡献值。首先，切换到教程目录：

1. #> cd 020-mmpbsa\Decompose

复制代码

残基能量分解的计算脚本和方法与上述示例完全相同。唯一区别在于残基能量分解和配体分子的精确选择，这需要修改以下参数：

1. -complex_selection         "protein and resid 165 or resname LIG" \  
   
2. -receptor_selection         "protein and resid 165"   \  
   
3. -ligand_selection           "resname LIG" \

复制代码

注："-complex_selection"、"-receptor_selection"和"-ligand_selection"参数用于选择复合物、受体和配体的分子部分，其功能类似于VMD Tk控制台中的"atomselect"命令。

运行以下命令可基于MM/PBSA快速计算残基自由能贡献：

1. #> molaical.exe -mmpbgbsa -f dec_mmpbsa_apbs.vmd

复制代码

输出结果包含每个残基的△G结合自由能，如下所示：

------------------------------------------------------------------------------------------------------------------------

   Sol:                  25                       0.4011                   1.9851                   9.9256

   Pol:                  25                       -0.5707                  1.9595                   9.7977

   Npol:                25                       -6.6169                  0.2676                   1.3378

   G binding:       25                       -7.1876        +/-  1.9791                   9.8957

------------------------------------------------------------------------------------------------------------------------

*所有能量值的单位均为千卡/摩尔 (kcal/mol)

注：

• 用户可采用上述类似方法，通过 MolAICal 计算残基对（pairwise）和单残基（per-residue）的自由能贡献。不同操作系统或 NAMD 版本可能导致结果存在差异，但若保留两位小数，结果将完全相同或高度相近。
• "输入文件与结果分析"的操作方式与上述结果高度相似，请参考前文说明。
  

总体而言，若参数预设得当，MolAICal 提供了一种快捷的 MM/PBSA 计算方法——仅需输入一条命令即可完成计算。

第四部分：熵的计算

3.3. 通过Carma和MolAICal进行熵的计算

3.3.1. 下载 Carma

Carma可以在Windows或Linux操作系统上运行。从https://github.com/glykos/carma, https://github.com/glykos/carma, http://utopia.duth.gr/~glykos/progs, 或 http://utopia.duth.gr/~glykos 上下载新版本Carma。

推荐使用最新版本的Carma进行熵的计算，之前的版本会出现“NaN”的错误等，具体请查看下面的链接：

https://groups.google.com/g/carma-molecular-dynamics/c/KpyY5sEkrj4

本教程选取100帧轨迹进行熵的计算。用户可以根据自己的实际情况选择合适的帧数，越多帧数越耗费计算时间。  

根据Carma的许可证（https://github.com/glykos/carma/blob/master/LICENSE），该许可证未限制使用、复制、修改、合并、发布、分发等权利。因此，MolAICal可直接集成并调用Carma来拟合分子并计算分子熵，而无需用户再次下载Carma。

选项：此处诸如'mpro.dcd'和'mpro.psf'等输入文件包含蛋白质、配体、离子、水分子等组分，可能占用大量计算机内存。熵值计算可能耗费较多时间和内存。为节省内存并提升计算速度，用户可通过以下命令分别提取复合物（受体与配体）、受体及配体的轨迹：

‌1. 复合物（第1分子）‌：

1. #> molaical.exe -call run -n vmd -c stripDCD -i -dispdev  text -psf "mpro.psf" -args protein,or,resname,LIG  "mpro.dcd" "complex" mpro.psf mpro.pdb

复制代码

注："complex"为输出文件前缀，程序将参照mpro.psf和mpro.pdb生成前缀为"complex"的文件，最终输出：'complex.psf', 'complex.pdb'和'complex.dcd'。

‌

2. 受体（第2分子）‌：

1. #> molaical.exe -call run -n vmd -c stripDCD -i -dispdev  text -psf "mpro.psf" -args protein "mpro.dcd"  "protein" mpro.psf mpro.pdb

复制代码

注："protein"为输出文件前缀，程序将参照原文件生成：'protein.psf', 'protein.pdb' 和 'protein.dcd'。

‌3. 配体（第3分子）‌：

1. #> molaical.exe -call run -n vmd -c stripdcd -i -dispdev  text -psf "mpro.psf" -args resname,LIG "mpro.dcd"  "ligand" mpro.psf mpro.pdb

复制代码

注："ligand"为输出文件前缀，生成文件包括：'ligand.psf', 'ligand.pdb', 和 'ligand.dcd'。

3.3.2. 计算复合物的熵

1. #> cd 020-mmpbsa\entropy

复制代码

切换到文件夹“com”。这一步是去掉复合物的旋转和平移：

1. #> molaical.exe -call run -centropy -i -v -fit -force -atmid ALLID -segid A -segid C complex.dcdcomplex.psf

复制代码

注：生成的文件"carma.fitted.dcd"采用固定命名，且该文件仅包含"-segid"参数指定的对应内容。用户可依据此提示优化调整命令参数以满足研究需求。

1. -call:   当值为"run"时，将调用外部程序
   
2. -c:      计算方式。当值为"entropy"时，将通过Carma计算熵值
   
3. -i:      在"-i"后输入Carma相关命令。该命令行参数"-i"必须放置在其他参数(如"-c"参数)之后  
   
4. -v:     详细信息输出  
   
5. -fit:    从DCD帧中移除全局旋转和平移  
   
6. -force:  指示Carma在熵值变为未定义前停止计算  
   
7. -atmid:  基于原子名称的原子选择  
   
8. -segid:  基于片段标识符的原子选择(仅限dcd-psf文件)

复制代码

这一步是计算熵值：

1. #> molaical.exe -call run -centropy -i -v -cov -eigen -mass -force -temp 310 -atmid ALLID -segid A -segid C-last 25 carma.fitted.dcd complex.psf >& com.log &

复制代码

3.3.3. 计算受体的熵

切换到文件夹“rec”。这一步是去掉受体的旋转和平移：

1. #> molaical.exe -call run -c entropy -i -v -fit -atmidALLID -segid A protein.dcd protein.psf

复制代码

这一步是计算熵值：

1. #> molaical.exe -call run -c entropy -i -v -cov -eigen-mass -force -temp 310 -atmid ALLID -segid A -last 25 carma.fitted.dcdprotein.psf >& rec.log &

复制代码

3.3.4. 计算配体的熵

切换到文件夹“lig”。这一步是去掉配体的旋转和平移：

1. #> molaical.exe -call run -c entropy -i -v -fit -force-atmid ALLID -segid C ligand.dcd ligand.psf

复制代码

这一步是计算熵值：

1. #> molaical.exe -call run -c entropy -i -v -cov -eigen-mass -force -temp 310 -atmid ALLID -segid C -last 25 carma.fitted.dcdligand.psf >& lig.log &

复制代码

3.3.5. 计算熵值

请检查log文件: “com.log”, “rec.log” 和 “lig.log”, 并在文件“com.log”,“rec.log” 和 “lig.log”中找到Andricioaei或Schlitter的熵值。例如：可以在“com.log”文件中找到Andricioaei 或 Schlitter的熵值 (见下图):

切换到上一层目录，基于日志文件运行MolAICal的操作方法，‌需使用"com.log"、"rec.log"和"lig.log"文件，执行步骤如下:

1. #> molaical.exe -entropy -if -c com/com.log -r rec/rec.log -l lig/lig.log -t 310.0

复制代码

这个结果显示:

Entropy Type:Andricioaei

The entropy (T *delta S) = -25.3034 (kcal/mol)

Entropy Type:Schlitter

The entropy (T *delta S) = -32.3195 (kcal/mol)

1. 各参数定义如下:
   
2.   
   
3. -entropy:   表示进行熵变计算（ΔS）  
   
4. -c:     复合物熵值或其结果文件（此处应为结果文件）  
   
5. -r:     受体（或第一分子）熵值或其结果文件（此处应为结果文件）  
   
6. -l:     配体（或第二分子）熵值或其结果文件（此处应为结果文件）  
   
7. -t:    分子动力学模拟的开尔文温度  
   
8. -if:   若输入参数包含"-if"，则从结果文件中读取熵值

复制代码

注意：Windows下计算熵值速度很慢，可能只有Schlitter的熵值；推荐使用linux系统进行熵的计算。推荐在Linux系统下计算熵值。

如果考虑熵值，MM/PBSA 通过以下方程计算：

1. ΔGbind = ΔH – TΔS ≈ ΔEMM+ ΔGsol - TΔS
   
2. ΔEMM=ΔEinternal + ΔEele +ΔEvdw
   
3. ΔGsol=ΔGPB + ΔGSA
   
4. ΔGSA= γ* SASA + β

复制代码

正如上文所述，熵能值未被纳入考虑，具体如下：

delta G binding: -66.5006     +/- 2.7810 (kcal/mol)

• 使用 Andricioaei 熵
  

ΔGbind = -66.5006 – (-25.3034) = -41.1972 (kcal/mol)

• 使用Schlitter熵
  

ΔGbind = -66.5006 – (-32.3195) = -34.1811 (kcal/mol)

注：若分子动力学（MD）模拟过程中未发生结合诱导的结构变化，或不同系统间的结构变化高度相似，则可省略熵计算‌，此时仅使用不考虑熵值的结合自由能（ΔG）数值就足够了。