自由能专题2：计算与分析指南

摘要：

基于分子动力学模拟的自由能计算显示了从药物发现到物理性质预测和结构功能研究的巨大应用前景。。 但是，这些计算仍然是困难和繁琐的分析，而且分析的最佳实践也没有得到很好的界定和推广。 基本上，分析这些计算的每个小组都需要决定如何进行分析，并且，通常，开发自己的分析工具。在这里，我们回顾并推荐从分子模拟中获得可靠自由能的最佳分析实践。 此外，我们提供了一个Python工具，alchemical-analysis.py，在GitHub上免费提供，作为pymbar包（位于http://github.com/choderalab/pymbar）的一部分（alchemical-analysis.py现已独立，在githup上：），它实现了这里所回顾的几个参考模拟包的分析实践，它可以适应处理其他包的数据。这篇综述和工具都涵盖了一般的化学计算分析，包括通过热力学积分和基于自由能扰动的估计器的自由能估计。我们的Python工具还可以处理多种类型的自由能计算的输出，包括扩展的集合和Hamiltonian复本交换，以及标准的fixed集合计算。我们还调查了一系列评估数据质量和自由能估计的统计和图形方法，并在我们的工具中提供了这些方法的原型。我们希望这个工具和讨论能成为自由能计算分析最佳实践的更多标准化和协议的基础。

1.引言：

1.1.自由能计算有助于药物研发

复杂的化学和生物系统对现代分子和计算科学提出了关键的挑战。我们寻求能够提供定量预测而不仅仅是定性见解的计算模型。研究人员试图回答诸如 “多少？”、“多大？”、"多紧？"等问题，并越来越多地应用物理细节计算来帮助回答这些问题。模型试图模仿或模拟有关过程，帮助揭示和提供对机制和现象的新理解，而这些机制和现象可能是具有挑战性的或不可能通过实验探究的[20, 23, 38, 43]。

自由能计算[7,9,11,18,25]是计算技术的一个很好的例子，它为一个特定的问题提供了定量的答案–在这种情况下，"‘系统的两个热力学最终状态之间的自由能差是多少？’"。例如，这个问题出现在药物发现[16]中，在这种情况下，需要根据药物与靶蛋白的结合力[9]等标准对药物进行排名。为了使用自由能计算来回答这个问题,我们必须首先建立一个系统模型,然后确定要计算自由能差的末端状态。

1.2.自由能的计算从结束状态的定义开始

以下并不需要模拟分子从气象中到水里整个过程，在气象中可以开关分子的相互作用即可起到相同的作用

热力学端态（例如图1中的状态1和2）是自由能计算的关键起点。原则上，端态之间的自由能差异可以简单地从一个或两个状态下进行的模拟计算出来[11]。 但在实践中，这对于生物分子系统来说，在合理的时间尺度上通常是不可能的。为了计算出准确的状态之间的自由能差，它们的相空间积分必须有足够的重叠，而这在实践中只有在两个状态极其相似的情况下才能实现。当不是这种情况时，就不可能直接计算出端态1和2之间的自由能差。在这种情况下，我们可以计算一系列中间态之间的自由能差，这些中间态有足够的重叠，从状态1通向状态2。这些中间状态通常是人为的，非物理状态，被构造成连接感兴趣的物理状态，并形成连接感兴趣的两个最终状态的热力学循环的一部分（见图1）。对于大多数与结合度、溶解度和溶解度相关的自由能计算，这种替代的、非物理的途径涉及有效地删除和/或插入一些原子，同时也可能对这些原子和其他原子进行参数改变。

​ （图1标准分子溶剂化自由能计算的热力学循环。）

这里，蓝色背景代表水，透明背景代表气相。

目标是求出两种状态之间的自由能差$\Delta G_{hydration}$：端态1（左上角）代表气相中的溶质，端态2（右上角）代表溶解的分子。

$\Delta G_{hydration}$是作为沿替代途径1引入的末端态1和2与中间化学态、中间态3(左下)和中间态4(右下)之间自由能变化的总和来计算的1-> 3 -->4 -->2，其中可能包括插值在1 -->3, 3 --> 4, 和 4 --> 2.

底部的溶质分子黑白相间的外观，表明溶质处于与环境没有非键相互作用的状态，也可能内部也没有非键相互作用。

溶质分子的黑白外观，在底部，表示溶质处于与环境没有非键性相互作用的状态，也可能没有内部非键性相互作用。这些情况分别称为解耦和湮灭，在文中讨论。在解耦的情况下，转化途径减少为单分支2 -->4，溶剂化自由能是这一路径自由能的负值。在这里所描述的湮灭的情况下，它也修改内部溶质相互作用，两条分支是必要的。

湮灭和解耦可以被认为主要是在如何处理电荷和LJ参数方面有所不同。具体来说，湮灭涉及实际设置溶质部分电荷为零，而解耦涉及关闭与环境的电荷相互作用。同样，湮灭涉及实际设置Lennard-Jones参数为零，而解耦涉及关闭与环境的相互作用。

我们目前的解释是专门针对更一般的情况，但在解耦的情况下，无需进行气体转化（1-> 3），并且总转化率降低为单个分支（2-> 4），水合自由能的变化被发现为$G_{decoupling}^{Water}$(与水解耦)的负值，但根据所采用的实验参考状态，可能存在理论标准状态修正的例外。

1.3.热力学循环描述了两种最终状态之间的交替路径

标准水化自由能计算的热力学循环是由连接四个感兴趣的状态的四条分支组成的（图1）。(2)一个分子与一盒水相互作用，(1)同一分子单独存在于气相中，(4)分子在水中但不与周围的水相互作用，(3)不相互作用的分子又单独存在于气相中。因此，我们通过修改溶质分子在每一个环境中的溶解自由能来计算。为此，我们计算关闭溶质与其环境的非键相互作用（称为解耦）或同时关闭内部非键相互作用和与环境的相互作用（称为湮灭）的自由能（图1）.1具体而言，我们计算与图1相关的自由能，1->3和2->4，分别关闭分子在气相和水中的相互作用。 最常见的是，这两种转换都是通过首先将溶质电荷缩放（通常是线性的）为零，然后关闭溶质的伦纳德-琼斯（LJ）相互作用（通常通过’‘软核’'方案[2]）来进行的。 这些转换是通过引入一个参数$\lambda$来调节系统的势能，所以当$\lambda$从0到1时，势能在初始状态和目标最终状态之间转换，通过一系列的步骤完成。然后，在连接这两个状态的一组不同的$\lambda$值上运行模拟。换句话说，两种转化途径中的每一种都被细分为各种单独的步骤，其中每一步都涉及两个$\lambda$值之间的转换。

选择$\lambda$值的数量和间距是为了保证被考虑的两种状态的构象空间之间有足够的重叠。这些中间状态及其相应的$\lambda$值通常被称为“炼金”，因为它们对应于非物理状态，往往涉及所考虑的物种的化学特性的改变。

在气相和溶液中分子发生化学转化后，仍需连接两个端态[图1中的（3）和（4）]。 但是，非相互作用分子的自由能并不取决于其环境的性质，所以非相互作用分子在环境之间的转移自由能为0。因此，从图1中的状态(3)到(4)没有相关的自由能变化。

1.4.相互作用是可以解耦的

静电相互作用的变化通常与LJ相互作用的变化分开，以避免在自由能量估计和采样挑战的不准确。具体来说，如果在原子的LJ相互作用被消除的同时，保留了静电相互作用，那么相关的电荷就会越来越多地暴露出来，并会产生巨大的静电力，导致巨大而昂贵的自由能差被收敛起来。在极端情况下，这可能导致正电荷和负电荷之间缺乏分离，这一个是特别的问题，可能导致数值不稳定和模拟崩溃[3,31]。分离这些变换的另一个好处是，它提供了一种机制来维持电荷相互作用[27]与k的最佳有效线性缩放，同时对LJ相互作用使用替代的缩放方案。具体来说，为了避免导数$\frac{(\partial U)}{(\partial \lambda )}$(TI分析所需dhdl.xvg文件)本来就不连续的情况，LJ变换的电势通常通过软核电势来处理[2, 36, 39, 42]。虽然我们在这里重点讨论的是解耦方案，但我们的一般分析也适用于组合的情况（不需要静电解耦的情况[5，12]）。

1.5.中间的 lambda 状态可以通过lambda向量来控制

GROMACS[33]和DESMOND[37]等MD软件包可以通过多个$\lambda$值控制不同相互作用类型的进度来处理自由能计算，这样就可以分别控制库仑、LJ和约束变换等自由能计算。沿着变换路径的每一步都与一组独特的$\lambda$值相关联，通常称为$\lambda$向量。对于图1中的热力学循环，k向量有两个组成部分，控制库仑和LJ相互作用。

每条变换路径1 -->3 和 2—>4可以呈现为一列中间耦合状态$({\lambda }_{coul}，{\lambda }_{LJ})$，初始和最终状态为(0，0)和(1，1)，如图2所示。如果$\lambda$控制了库仑和LJ与溶质环境相互作用的强度，那么随着$\lambda$进展，溶质-溶剂的相互作用逐渐减少，直到在最终状态下，系统由纯溶剂组成，被一个由无相互作用的、具有完全内部相互作用的孤立溶质组成的平行系统所覆盖。在计算中，内部溶质非键相互作用也被去除，溶质端态略有不同，由只通过其键相互作用的原子组合而成。

一旦选择了$\lambda$个状态，就进行平衡模拟，存储必要的分析信息（尤其是$\frac{(\partial U)}{(\partial \lambda )}$，电位相对于$\lambda$的导数，$\Delta U_{ij}$，不同$ \lambda$值下状态间的电位能差，来自别轨迹）

（图2：变换图1中的1–>3路径和2—>4路径）

该过程可以看作是$\lambda$向量空间从（0,0）(填充圆)到（1，1）(空心圆的)的一条路径。在笛卡尔坐标平面上，这条路径被显示为由横轴${\lambda }_{vdWaals}$和纵轴${\lambda }_{Coulomb}$形成的淡蓝色方块，这两条轴分别控制了范德华和库伦相互作用。

出于文中所讨论的原因，我们从以下几点入手 (0,0)，对应于完全相互作用的分子，而经过几个炼金中间状态(位置为用箭头表示）沿${\lambda }_{Coulomb}$轴（红色箭头）直到我们到达(0,1)点，这对应于静电的非相互作用的分子。然后，我们修改${\lambda }_{vdWaals}$（绿色箭头）直到我们达到目标状态(1,1)，对应于非相互作用的分子。如果换成控制转化进度而不是一个向量，变换路径为会沿着正方形的对角线(蓝色箭头)。

1.6自动分析应该是自由能计算的一个重要部分

自由能量计算通常是专家们的工作，原因之一是其设置和分析通常需要大量的人工干预。随着越来越多的研究人员参与到自由能计算中来，标准的分析工具变得越来越重要，既可以帮助确保最佳实践得到遵循，又可以避免重复劳动。

我们这里的重点是对自由能计算的分析，自由能计算通常包括一系列连续的步骤。这些自由能计算本身可以用各种不同的取样技术来进行，我们这里主要关注分析阶段，无论
采样技术的。

在这里，我们提出我们认为是目前分析炼金自由能计算的最佳实践。在概念上，我们将分析分为四个主要阶段。

1.对数据进行子抽样，保留不相关的样本。

2.通过各种基于TI和FEP的方法计算自由能差以及相应的统计误差。

3.对计算出的数据进行文字和图形输出

4.检查是否有：

收敛性，并确定模拟的平衡部分。

对所有相邻的lambda状态对有良好的相位空间重叠。

2.概念、理论和自由能估算

我们将注意力集中在对模型自由能计算的分析上–在本例中，我们选择3-甲基吲哚的水合自由能计算作为说明性的例子。

后面的讨论将假设读者已经进行了自由能计算，并希望对结果数据进行分析。

在我们的案例中，我们已经在GROMACS中运行这个自由能计算，我们提供了一个Python工具，它实现了这里描述的GROMACS、SIRE（http://siremol.org/Sire/Authors.html）和AMBER[6]数据文件的程序，并提供了例子。

然而，除了读取输入数据外，我们的代码与具体的模拟包无关，并且可以很容易地适应于包含这里使用的自由能量估计器所需数量的任何数据格式（$\frac{(\partial U)}{(\partial \lambda )}$用于TI，和$\Delta U_{i.j}$用于FEP计算)。因此，虽然我们这里的例子使用的是GROMACS仿真包的情况，但我们的原型工具，可以在GitHub上免费获得https://github.com/choderalab/pymbar-examples，可以很容易地修改，以便与其他仿真包一起工作。

2.1.获取输入数据

如前所述，为了全面的计算自由能、一般分析所需的关键输入信息包括（至少）相邻的$\lambda$态之间的势能差，以及所有$\lambda$态的$\frac{(\partial U)}{(\partial \lambda )}$值。具体来说，我们将给出一个使用GROMACS格式的输入文件进行计算的例子。特别是，GROMACS 目前（从 v3.3 到 v5.0）将所有能量存储到二进制能量文件中，但也将分析所需的所有势能及其差值写入人可读的 .xvg 格式的文本文件中。

在GROMACS中，这些文件的格式如图3所示。对于标准模拟(与下面讨论的扩大合奏模拟相反)，有几个这样的文件–每$\lambda$值都有。在GROMACS中，$\frac{(\partial U)}{(\partial \lambda )}$字段的精确数量随着被利用的不同类型的k值的数量而变化，这与$\lambda$向量的维数相对应。

(图3. 3-甲基吲哚在水中的前十个快照的GROMACS自由能计算数据样本，数据出现在dhdl:1:xvg文件中，标题中给出了一行中每个字段的名称。这些都是：时间（皮秒），系统的能量（根据所使用的选项，可以是势能或总能量），相对于扰动（$\frac{\partial U}{\partial {\lambda }_{Coul}}$和$\frac{\partial U}{\partial {\lambda }_{vdW}}$）中所有$\lambda$类型的总能量导数，评价指数1的当前$\lambda$状态和其他状态($\Delta U_{i=1,j,j=0,1,2...10}$)之间的总能量差，当没有质量扰动时，它被还原为势能的差异。从GROMACS 4.0版本开始，默认设置是只评估相邻状态($\Delta U_{1,0},\Delta U_{1,1},\Delta U_{1,2}$)之间的$\Delta U_{i,j}$，每当需要评估所有差异时(mbar计算方法)，应该设置.mdp文件中的 calc-lambda-neighbors = -1 （计算每一个状态的左边和右边的状态）。)

扩大的集合模拟[21,22,24,26,29]是一种方法，它允许同时探索k和坐标空间在一个单一的模拟，潜在地允许更快的采样跨化学状态，前提是划分构象相关的交替$\lambda$状态的动力学障碍不是很大。在扩展的集合模拟中，单次模拟对所有状态进行采样，从而产生一个单一的能量文件。

在每个时间步骤，模拟是在一个特定的$\lambda$状态，仅在扩展的合奏情况下，该状态才存储在输出能量文件的第二个字段中。这样就可以确定存储的$\frac{(\partial U)}{(\partial \lambda )}$和$\Delta U_{ij}$值属于哪个$\lambda$状态。

2.2.炼金分析技术可分为以下几个系列

如前所述，有多种方法可以从炼金自由能计算中获取输出，并产生自由能差。从概念上讲，这些方法可以根据用于计算$\Delta G$的数量分为两类：热力学积分(TI)17方法和自由能扰动(FEP)44方法。在TI中，沿$\lambda$状态组成的路径上的自由能变化被计算为势能函数导数与耦合参数$\lambda$的集合平均值的加权和：

$$\Delta G=\sum _{i=1}^K{W}_i⟨\frac{\partial U}{\partial \lambda }{⟩}_{{\lambda }_i}\text{\hspace{1em}(1)}$$

其中$W_i$为加权系数，取决于所使用的数值积分方案[28]。

TI中的数值积分有几种不同的方案。在我们提供的工具alchemical-analysis.py中，我们实现了TI-1和TI-3[28]，它们在数据点之间的插值方式上有所不同。TI-1使用梯形规则（一阶多项式），而TI-3使用（自然）立方花键。这些不同的TI方法的相对性能将取决于基础数据的性质和被积分的$\frac{\partial U}{\partial \lambda }$曲线的形状–因此，它取决于炼金术的选择的路径。

基于扰动的方法包括一系列与FEP松散相关的技术。在我们的原型工具中，这些技术包括删除指数平均法（Deletion Exponential Averaging,DEXP）、插入指数平均法（Insertion Exponential Averaging,IEXP）、高斯删除法（Gaussian Deletion,GDEL）、高斯插入法（Gaussian Insertion,GINS）、Bennett接受率（Bennett Acceptance Ratio,BAR）、非优化Bennett接受率（Unoptimized Bennett AcceptanceRatio,UBAR）、基于范围的Bennett接受率(Range-based Bennett Acceptance Ratio,RBAR)和多态Bennett接受率(Multistate Bennett Acceptance Ratio，MBAR)。其中一些，如BAR[1]和MBAR[35]，已被普遍使用，并在本领域的术语中根深蒂固，而另一些则没有习惯的、被普遍接受的名称，或者仍然鲜为人知。对于这些方法我们将使用建议的命名惯例 Paliwal和Shirts的报告[28]。

前两种，DEXP 和 IEXP，是基于两个相邻状态之间势能的指数平均方案–所谓的 Zwanzig 关系[44]：
$$\Delta G_{ij}=-\frac{1}{\beta }ln⟨exp(-\beta {\Delta }_{ij)}{⟩}_i\text{\hspace{1em}(2)}$$
根据转换的方向，该过程可以被解释为’‘删除’‘或’‘插入’’，因此缩写中的第一个字母。通常情况下，其中一个过程，DEXP，是朝着增加熵的方向进行的，而另一个过程，IEXP，是朝着减少熵的方向进行的。

如果势能差以高斯方式分布(GINS和GDEL[28])，那么Zwanzig关系就简化为[13，14]:
$$\Delta G_{ij}=⟨\Delta U_{ij}⟩-\frac{\beta }{2}{\sigma }_{\Delta U_{ij}}$$