基于百度AI Studio的药物设计——NAMD自由能微扰计算（FEP笔记本）

药物设计中最重要的任务之一是在一系列先导候选药物中预测哪些药物与治疗靶点的结合更强。在这个方向上，已经开发了相对结合自由能方法，它依靠基于物理的分子模拟和严格的统计力学来计算母体候选药物和类似物之间结合自由能的差异。例如，自由能微扰 (FEP) 计算将初始（参考）分子和最终（目标）分子之间的自由能差计算为通过初始状态采样评估的能量差函数的平均值。

然而目前支持FEP的商业软件价格昂贵，学术软件有众多限制。NAMD作为一款强大且开源的动力学模拟软件，支持并行运算、GPU加速、多种能量计算方式、变量调节，可快速模拟大分子体系的分子动力学代码。结合百度提供的免费算力，或许每一个新药研究者都可以对自己的体系进行FEP计算，为设计提供思路，加速开发，节约成本。

项目介绍

  1. 项目背景

当与T4 Lysozyme结合的小分子甲苯变成苯时，其结合自由能发生了怎样的变化？

在这里插入图片描述

为预测计算上述问题，项目以

1）伊利诺伊大学开发的动力学模拟软件NAMD

http://zarbi.chem.yale.edu/ligpargen/

2）耶鲁大学开发的配体力场参数生成网站

https://feprepare.vi-seem.eu/indexlpg.php

3）雅典大学开发的基于NAMD FEP的准备工具FEPrepare网站

https://www.ks.uiuc.edu/Research/namd/

为基础（请酌情引用。这些网站上几乎包含了有关于动力学模拟的一切，包括FEP的理论知识），提供一个易入手的FEP计算流程。

使用者可以在准备好输入文件后，直接使用百度AI Studio的脚本任务，一键开启FEP之旅 https://aistudio.baidu.com/aistudio/projectdetail/4326308?contributionType=1 ；也可以在该笔记上，从头准备输入文件，了解NAMD FEP的方法https://aistudio.baidu.com/aistudio/projectdetail/4326115?contributionType=1 。

此外，本项目源于https://github.com/quantaosun/NAMD-FEP

在这里插入图片描述

感谢孙全涛博士的指导

  1. 数据集介绍

数据集1来自NAMD网站，NAMD_3.0alpha8_Linux-x86_64-multicore-CUDA.tar.gz（FEP执行者）

https://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=NAMD

数据集2来自FEPrepare网站，下载你的结果（FEP执行对象）

https://www.ks.uiuc.edu/Research/namd/

  1. 应用场景

基于结构的药物设计，比较一对小分子（binding pose相似）的结合自由能差异，即在药物设计中，从化合物A到化合物B，其活性好了多少，或差了多少。

注意：需要明确的binding pose！！

在这里插入图片描述

关于FEP的理论知识：计算ddG，可以直接模拟两次化合物从溶剂状态到嵌入蛋白口袋，再相减（即计算dGToluene-dGToluene）。但是，模拟化合物从溶剂状态到嵌入蛋白口袋是一个极其复杂且漫长的过程，耗时巨大且不准确。**通过如图的循环，我们可以将这个复杂过程转换为对计算来说较为容易的过程（dGcomplex-dGsolvent），因为热力学只在意其最终状态，不在乎过程中发生了什么。**这个过程之所以比较容易是因为它的状态变化小（从溶剂到溶剂，或从complex到complex），在模拟时更容易达到平衡，更省时，也更准确。所以在结果分析时，将计算所得的dGcomplex和dGsolvent相减，就能得到我们想要的ddG。

此外，FEP计算还涉及炼金术化学的概念，从甲苯到苯，我们不必考虑它是怎样转化的，不必在意在化学上的反应可行性；也不是将整个分子替换。只需要考虑发生变化的位置，模拟中逐渐关闭原分子变化的位置的作用力，打开新分子在该处的作用力。

这些在namd的User’s Guide中都有介绍，但namd目前只给出了很简单的教程，甚至没有涉及ligand-protein的体系，且从头学习namd的将耗费大量的精力，这也是本项目的诞生的起因。

方案流程

在这里插入图片描述

首先，你需要一个初始的complex.pdb文件，它可以来自PDB数据库，或者十分可靠的pose结果。

对complex.pdb中的配体做修改，设计你的化合物，并和原配体叠合对齐（align）。

分离蛋白于与配体，分别保存成3个pdb文件：蛋白用pdbfixer处理；配体用ligpargen处理。

处理完的文件用feprare得到FEP的输入文件。

上传输入文件后，在AI Studio中依次运行bash脚本。

运行结束后，下载四个fepout文件，用本地的VMD软件进行结果分析。

安装

!tar -xvf /home/aistudio/data/data119630/NAMD_3.0alpha8_Linux-x86_64-multicore-CUDA.tar.gz -C /home/aistudio
!yes |conda install -c conda-forge pdbfixer

Fetching package metadata .............
Solving package specifications: .

# All requested packages already installed.
# packages in environment at /opt/conda:
#
pdbfixer                  1.8.1              pyh6c4a22f_0    conda-forge
yes: 标准输出: 断开的管道

指定工作路径

#your work directory
%mkdir /home/aistudio/work/T4L
%cd /home/aistudio/work/T4L


mkdir: 无法创建目录"/home/aistudio/work/T4L": 文件已存在
/home/aistudio/work/T4L

0_输入文件

上传蛋白pdb文件（4w53.pdb），用pdbfixer预处理蛋白，避免FEPrepare失败

此步获得的protein_standard.pdb，与ligpargen获得的6个配体相关文件（共7个文件）上传至Fepreare，即可获得FEP输入文件

!pdbfixer 4w53.pdb --pdbid=4w53 --output=protein_standard.pdb --keep-heterogens=none --add-residues --replace-nonstandard --add-atoms=all

Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.

1_路径

执行此步前，确保你已经获得来自Fepreare的FEP输入文件，并上传至工作路径

在这个例子的中，我的FEP输入文件为ligand0_ligand1_07112022051655

!unzip ligand0_ligand1_07112022051655.zip #your FEP input file.
!bash /home/aistudio/1_pre.sh

Archive:  ligand0_ligand1_07112022051655.zip
  inflating: reference.rtf           
  inflating: top_opls_aam.inp        
  inflating: mutant.rtf              
  inflating: ligand.pdb              
  inflating: complex.pdb             
  inflating: par_all22_prot_metals.inp  
  inflating: mutant.pdb              
  inflating: updated.prm             
  inflating: complex/complex_wb.pdb  
  inflating: complex/top_opls_aam.inp  
  inflating: complex/vmd_prepare_complex_after_gui_autopsf_ionize  
  inflating: complex/complex_wb.psf  
  inflating: complex/complex-input-files.tar.gz  
  inflating: complex/ligand.pdb      
  inflating: complex/psf-complex.pdb  
  inflating: complex/ionized.fep     
  inflating: complex/min-max_center  
  inflating: complex/complex.pdb     
  inflating: complex/psfgen          
  inflating: complex/ionized.pdb     
  inflating: complex/vmd_prepare_complex_after_gui_autopsf  
  inflating: complex/ionized_complex.fep  
  inflating: complex/ionized_complex.pdb  
  inflating: complex/complex_wb.log  
  inflating: complex/vmd_log.txt     
  inflating: complex/ionized.psf     
  inflating: complex/chainX.pdb      
  inflating: complex/psf-complex.psf  
  inflating: complex/chainA.pdb      
  inflating: complex/ligand.rtf      
  inflating: solvent/top_opls_aam.inp  
  inflating: solvent/ligand_wb.pdb   
  inflating: solvent/psf-solvated.psf  
  inflating: solvent/ligand.pdb      
  inflating: solvent/ionized.fep     
  inflating: solvent/min-max_center  
  inflating: solvent/complex.pdb     
  inflating: solvent/ligand_wb.log   
  inflating: solvent/ligand_wb.psf   
  inflating: solvent/vmd_prepare_ligand_after_gui_autopsf  
  inflating: solvent/psfgen          
  inflating: solvent/vmd_prepare_ligand_after_gui_autopsf_ionize  
  inflating: solvent/ionized.pdb     
  inflating: solvent/psfgen_solv     
  inflating: solvent/ionized_solvent.pdb  
  inflating: solvent/psf-solvated.pdb  
  inflating: solvent/ionized_solvent.fep  
  inflating: solvent/solvent-input-files.tar.gz  
  inflating: solvent/vmd_log.txt     
  inflating: solvent/ionized.psf     
  inflating: solvent/ligand.rtf      
  inflating: README-FEPrepare.pdf    
  inflating: fep.tcl                 
  inflating: top_all22_prot_metals.inp  
  inflating: reference.pdb           
  inflating: mutant.prm              
  inflating: reference.prm           
  inflating: ligand.rtf              
  inflating: par_opls_aam.inp        
md_backward_10.namd
md_backward_11.namd
md_backward_12.namd
md_backward_13.namd
md_backward_14.namd
md_backward_15.namd
md_backward_16.namd
md_backward_1.namd
md_backward_2.namd
md_backward_3.namd
md_backward_4.namd
md_backward_5.namd
md_backward_6.namd
md_backward_7.namd
md_backward_8.namd
md_backward_9.namd
md_forward_10.namd
md_forward_11.namd
md_forward_12.namd
md_forward_13.namd
md_forward_14.namd
md_forward_15.namd
md_forward_16.namd
md_forward_1.namd
md_forward_2.namd
md_forward_3.namd
md_forward_4.namd
md_forward_5.namd
md_forward_6.namd
md_forward_7.namd
md_forward_8.namd
md_forward_9.namd
npt_equil.namd
nvt_equil.namd
md_backward_10.namd
md_backward_11.namd
md_backward_12.namd
md_backward_13.namd
md_backward_14.namd
md_backward_15.namd
md_backward_16.namd
md_backward_1.namd
md_backward_2.namd
md_backward_3.namd
md_backward_4.namd
md_backward_5.namd
md_backward_6.namd
md_backward_7.namd
md_backward_8.namd
md_backward_9.namd
md_forward_10.namd
md_forward_11.namd
md_forward_12.namd
md_forward_13.namd
md_forward_14.namd
md_forward_15.namd
md_forward_16.namd
md_forward_1.namd
md_forward_2.namd
md_forward_3.namd
md_forward_4.namd
md_forward_5.namd
md_forward_6.namd
md_forward_7.namd
md_forward_8.namd
md_forward_9.namd
npt_equil.namd
nvt_equil.namd
nvt-equilibration/nvt_equil.namd
npt-equilibration/npt_equil.namd

2_PBC

自动填写周期性边界条件

!bash /home/aistudio/2_PBC_build.sh

the x-axis of your system is：
cellbasisvector1         93     0     0
the y-axis of your system is：
cellbasisvector2         0     71     0
the z-axis of your system is：
cellbasisvector3         0     0     83
the center of your system is：
cellOrigin          -14.8 14.3 12.3
the x-axis of your system is：
cellbasisvector1         34     0     0
the y-axis of your system is：
cellbasisvector2         0     35     0
the z-axis of your system is：
cellbasisvector3         0     0     36
the center of your system is：
cellOrigin          -33.3 6.2 4.0

3_配置

在这个例子中，只修改了模拟时长（nvt\npt都为0.3ns，每个λ为1ns，即总时长65.2ns；在AI Studio的脚本任务中，选择P40环境，约耗时12个小时）；

当然如果需要，可以进行更多细节上的修改

!bash /home/aistudio/3_runtime_setting.sh

4_运行

!bash /home/aistudio/4_run.sh >> /home/aistudio/namd_fep.log

^C

5_结果处理

合并1至16.fepout文件

!bash /home/aistudio/5_merge.sh

6_结果分析

http://www.ks.uiuc.edu/Research/vmd/plugins/parsefep/

所得的四个complex_forward.fepout、complex_backward.fepout、solvent_forward.fepout、complex_backward.fepout，借助vmd（本地）的ParseFEP插件进行分析，本例结果如下

分析前需要了解一些FEP理论知识，以本次计算为例，FEP将甲苯到苯的过程称为划分为16个λ，每个λ为0.0625个过程，每0.0625个过程模拟时长为1ns；且整个过程需要在complex和solvent两种环境中各运行两次（正着跑一次，“甲苯到苯”，逆着跑一次“苯到甲苯”）；即共64ns。所有的结果整合一起，误差修正后，将complex的值减去solvent的值，得到最后结果ddG。（如果对这些分析不感兴趣，直接看ddG结果就好~）

complex:

在这里插入图片描述

此图是对complex部分的FEP总览，可以在0.5个过程后dG的变化平稳了一些。黑线为正着跑，红色为逆着跑，两者较好地重合。继续延长模拟时间，可增加它的稳定程度。但要得到最后精确度，往往需要大量的时间，对于本例子这个较简单的体系，这个程度已经够了。
在这里插入图片描述

此图为complex部分每个λ的dG展示，黑线为正着跑，红色为逆着跑，两头平稳汇与一条平线时为好。不好的λ意味着对计算误差有较大的影响，该λ需要修改，或延长时间。

在这里插入图片描述

此图为complex部分每个λ的dG概率分布，即误差范围，同样黑线为正着跑，红色为逆着跑。

SOS-estimator: total free energy change is 9.876073900572676 , total error is 0.04620452645958048

complex部分的最终计算结果。

solvent：

解释可参考complex部分
在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

SOS-estimator: total free energy change is 9.596021229259472 , total error is 0.013017820673635166

ddG

= 9.876073900572676 ( ± 0.04620452645958048 ) - 9.596021229259472 ( ± 0.013017820673635166 )

= 0.280052671 ( ± 0.059222347 ) kcal/mol

甲苯到苯的结合自由能之差为 0.28 kcal/mol，即计算结果为 苯结合T4L的活性比甲苯差 0.28 kcal/mol。

相比于实验值，十分接近！

ddG实验值

= 0.33 kcal/mol

在这里插入图片描述

在这里插入图片描述

实验值，Figure 8，以及Table 4 来自参考文献：

Pall, M. (2013). CHARMM-GUI Free Energy Calculator for Absolute and Relative Ligand Solvation and Binding Free Energy Simulations. International Journal of Molecular Sciences, 14(11), 22274–22330. https://doi.org/10.1021/acs.jctc.0c00884.CHARMM-GUI

此文仅为搬运，原作链接：https://aistudio.baidu.com/aistudio/projectdetail/4326115?channelType=0&channel=0

# 百度 # 人工智能 # php
Logo
百度飞桨AI Studio社区

学大模型，用大模型上飞桨星河社区！每天8点V100G算力免费领！免费领取ERNIE 4.0 100w Token >>>

更多推荐

cover

利用Amazon Bedrock畅玩Claude 3等多种领先模型，抢占AI高地(体验倒计时4小时)

avatar 百度飞桨星河社区
cover

RAPTOR：索引树状 RAG，使用树结构来捕捉文本的高级和低级细节

avatar 百度飞桨星河社区

RAG介绍

LlamaIndex是一个基于大型语言模型（LLM）的应用程序数据框架，旨在连接LLM与外部数据，特别是私有或特定领域的数据，以实现更准确的文本生成和知识增强。

avatar 百度飞桨星河社区