目的：采用Discovery Studio，以一组小分子化合物为实例，示范FEP计算及结果分析操作过程。

所需功能：Discovery Studio client，DS Preparing analogs and ligand pairs、DS Setting up the relative FEP calculations、DS Running the FEP calculations、DS Collating and analyzing the FEP results

所需文件：4YKR_prep.dsv、hsp90_analogs.sd、hsp90_lead.sd

所需时间：30分钟+FEP计算时间（与电脑硬件相关）

介绍：理解和量化配体和蛋白质之间相互作用强度的能力在药物发现项目中是至关重要的，以识别配体的变化如何可能改变亲和力。对结合亲和度的可靠估计可以改进先导优化过程中配体的选择。在早期先导优化阶段，可能会考虑成千上万种不同支架家族的配体，因此诸如对接和MM-GBSA等计算方法适用于所需的吞吐量。后期先导优化涉及小的、顺序的化学修饰，需要更精确，可能需要明确的溶剂建模和严格的自由能方法。炼金术自由能摄动(FEP)是建立在热力学循环基础上的。

水平物理过程描述了蛋白质复合物中从溶剂中未结合的配体到结合态的自由能变化。在垂直炼金术过程中，配体1在非结合态和结合态转变为配体2。这是通过一系列非物理中间状态进行计算的。ΔG_complex - ΔG_free的炼金术相对束缚自由能的计算相当于ΔG2 - ΔG1，但计算成本较低。

本教程主要包括以下步骤：

• 制备类似物和配体对

• 建立相对FEP计算

• 运行FEP计算

• 整理和分析FEP结果

一、制备类似物和配体对

FEP动力学模拟需要所有感兴趣的配体定位在蛋白质活性位点上。在第一步中，您将使用生成模拟构象协议将您的配体与蛋白质活性位点中的先导配体对齐。该协议通过最大公共子结构(MCS)方法检测每个模拟物和先导化合物之间的公共原子和不同原子。它直接从相应的铅原子分配模拟物的公共原子的坐标，并为不同的原子生成新的坐标。

热休克蛋白90蛋白配体复合体的PDB编码为4YKR。本教程中使用的示例文件是从PDB网站下载的，蛋白质是用Prepare Proteins协议制备的。在制备蛋白质时，要特别注意晶体水，以确定它们是否应该作为保守的水分子保留。深埋在蛋白质核心的水分子可能在维持蛋白质的折叠结构中发挥作用，并经常与极性基团形成强氢键。在本教程中，我们在结合位点保留了三个已知的与先导化合物形成重要相互作用的水分子，并在整个工作流程中使用相同的蛋白质结构。然而，在未来的应用中，特定的水分子的相关性最初可能不知道，建议生成模拟构象协议运行时，蛋白质中不存在任何晶体水分子，FEP设置与蛋白质结构保留所有的水分子。

1. 打开4YKR_prep.dsv、hsp90_analogs.sd、hsp90_lead.sd这三个文件，位于Samples | Tutorials | Simulation |文件夹。

2. Generate Analog Conformations，该程序根据一系列先导化合物的一个或多个先导化合物的三维几何结构和相应的目标蛋白质结构，允许获得同型模拟化合物系列中配体的三维几何结构。该协议首先通过最大公共子结构(MCS)方法检测每个模拟物和先导化合物之间的公共原子和不同原子。它直接从相应的铅原子分配模拟物的公共原子的坐标，并生成不同原子的坐标。由于模拟化合物可以有不止一种方式与先导化合物拓扑匹配，而且由于不同的原子会生成多个构象，因此输出通常会导致每个模拟化合物产生多个构象。如果提供了蛋白质受体，对于每一个构象，通过对复合体和自由蛋白质和配体进行单点CHARMm能量计算得到结合能，并以能量最有利的构象优先排序。然而，这种能量不太可能可靠地预测结合亲和性，因为它不涉及结构优化。打开Generate Analog Conformations按照下图进行参数设置，运行。

3.  结果文件中Analog Conformations拥有8个构象，其中hsp90_52有两个构象，Top Scoring Conformations文件自动集合和打分最高的每个配体的构象，因此选择这个文件即可。给该7个小分子添加力场。Simulation | Change Forcefield，选择CHARMm力场，Apply Forcefield应用力场（FEP算法同时支持CHARMm和charmm36两种力场，若要使用charmm36则用Type Ligands with MATCH (Prototype)程序应用）。

4. 打开Generate Ligand Pairs for FEP....，该程序的工作原理是首先为每个唯一的配对计算一个包含最大公共子结构(MCS)的相似性矩阵。从这个矩阵中选择配体之间的连接，使大多数相似的配体连接到可能的最佳范围，这是使用最大相似度算法完成的。对于含有N个配体的数据集，结果是N-1对的集合。 按照如下设置，因为第3步的hsp_90 lead文件就是hsp_12，因此将在Top Scoring Conformations文件里的hsp_12选择为先导配体（可以自行指定），运行。

5. 运行结果如下图所示，

分子窗口中的配体对被定义为第一个配体通过炼金术转变为第二个配体的组合。网络图显示了配体对之间的关系，其中节点是配体，边表示配对。边缘被标记为一对的相似度分数，箭头表示炼金术变化的方向。你可以选择指定一个先导配体，在这种情况下，每对配体根据它们与先导配体的距离排列，其中距离是由配体和先导配体之间的介入节点的数量来衡量的。hsp90_12配体与在4YKR中结晶的配体相同。PDB复杂，有了网络图，你可以从先导开始，向外单向移动到越来越不相似的类似物。

6.  接着我们需要准备运行FEP计算的结构。根据在拓扑和非键参数中识别出的共有和不同原子集，每个配体对生成合并的对偶拓扑分子。还生成了这种合并的对偶拓扑分子的溶剂化配体体系和溶剂化蛋白质配体络合物。打开Set Up Relative FEP Calculations…按照如下进行设置，我们可以设置指定到周期性边界条件的最大距离以及cell shape形状、是否添加盐溶液环境等，运行。

7.  在炼金术FEP计算中，利用合并的对偶拓扑分子将一个配体转化为第二个配体。这种嵌合分子是通过最大共子结构(MCS)方法首先检测配体之间的共有原子和不同原子而产生的。普通原子是那些在拓扑结构、原子类型和部分电荷方面被认为是相同的原子。所有其他原子都被指定为不同的原子组。运行结果如下，其中一个包含合并的对偶拓扑分子在一个立方形状的盒子里，水分子和盐离子。另一个分子窗口包含了每个对偶拓扑分子与蛋白质在类似溶剂化盒子中的复合物。

8. 分界点。前期的文件准备工作到此结束。正式的FEP计算在DS中有CPU和GPU（截止版本为DS2022，暂时只支持Linux系统）两种方式。下面分别展示。

9.1.  CPU计算。运行两遍Free Energy Perturbation (FEP)程序，一遍为计算溶剂化受体，一遍计算溶剂化配体，如下图所示。我们可以在这里去修改一些参数诸如Lambda窗口和生产的时间步长等以适应不同的计算需求。

运行结果如下：

其中一个FEP结果报告。

10.1 整理和分析FEP结果，在资源管理器中打开上述结果文件（也可以在job中右键locate定位），其中包含一些文件，如下。

这些文件夹中的结果包含每个配体对的一个simulationinfo.xml文件和一个自由能数据文件(.fepout)。xml文件指定了每个扰动中涉及的配体的名称，以及模拟是涉及自由配体还是复配体。.fepout文件是在NAMD FEP生产运行期间生成的，其中包含原始自由能数据。

打开Collate FEP Simulation Results....程序，将FEP数据目录配体参数设置为溶剂配体上NAMD FEP计算的结果位置，将FEP数据目录复合体参数设置为溶剂复合体上NAMD FEP计算的结果位置。（如果需要计算绝对自由能，则展开绝对自由能参数组，并在先导配体名称参数中输入hsp90_12。在先导配体自由能参数字段中输入值-10.27（实验值）。设置输入配体参数为Top_Scoring_Analog_Conformation:All.），运行。

运行结果如下所示，列出一个计算出的配体对的相对自由能值的表，以及一个配体的绝对自由能表。后一个表只有在你提供了先导配体和自由能值时才会生成。这两个表也都用柱状图图形化表示。

你会注意到，hsp90_52的绝对自由能在表格中被标记为“可能不可靠”，在图表中被标记为红色。这是因为经验观察表明，在同属数列中，相对自由能值的大小很少大于3千卡/摩尔。极值阈值参数允许您标记这些可能不可靠的预测，这可能是由于模拟中的不准确性，如由于采样不足或力场不收敛。对于某些包含非活性分子的数据集，3千卡/摩尔的阈值可能标志着大量的预测。在这种情况下，通过将阈值提高到4或5千卡/摩尔来测试结果对阈值的敏感性是有用的。

如果你将这些值与报告中计算出的绝对ΔG值进行比较，你会发现，除了不可靠的hsp90_52预测外，自由能摄动计算能够合理地预测绝对值。

也许比重现精确的绝对自由能预测更重要的是，FEP方法能够对配体进行排序。当绝对自由能预测不能计算出排序时，你仍然可以使用相对值进行分析。看看配体对hsp90_40 - hsp90_4的ΔΔG值为2.864，看看hsp90_4 - hsp90_3的值为-0.482。你可以据此推断，hsp90_40的结合能比hsp90_4好，而hsp90_3的结合能比hsp90_4好。配体对hsp90_4 - hsp90_5的ΔΔG值为2.746，说明hsp90_5的结合能最差。

Collate FEP模拟结果协议还生成一个网络图(如图)，图形化地表示配体对的ΔΔG值。边缘用配体对的相对自由能进行标记和颜色编码，箭头表示炼金术变化的方向。

在本教程中，您已经使用了相对炼金术自由能扰动来为后期先导优化的一系列同型配体排序。模拟所需的计算成本限制了该方法的小数据集，但其提供更精确的结合能定量估计的潜力使其成为基于合理结构的药物设计的宝贵工具。

9.2 GPU计算。GPU计算目前仅能在Linux系统（DS服务器安装在Linux系统，DS客户端依旧可以安装在windows上）上运行，windows系统暂时不支持。GPU计算整合了分析结果，可以一步生成结果。打开CHARMm Relative FEP Calculations (GPU)程序，参数如下。GPU程序下，我们仅需跑一次程序即可，在这里可以如CPU程序一样修改一些参数诸如Lambda窗口和生产的时间步长等以适应不同的计算需求。同时如果需要计算绝对自由能，则选择先导配体，填写实验值即可。

结果如下：

值得注意的是，GPU运行FEP程序有三种FEP Systems可供输入，分别是Binding（需要溶剂化的复合物和溶剂化的配体作为输入）、Hydration（需要溶剂化配体和真空配体作为输入，可以结算水合自由能）、Single（运行三个可能的系统中的一个，在运行两个单一类型的计算后，手动结合结果，以获得Binding或Hydration自由能）。其中真空配体的文件可在Set Up Relative FEP Calculations的输出结果中获得。