结合自由能计算（四）
基于MM/PBSA自由能算法的操作实例
前面提到的基于MM/PBSA来计算结合自由能的方法是目前应用较为广泛的一种自由能计算方法，它属于基于主方程算法中一种比较成熟算法。
图1 MM/PBSA自由能计算的流程图
基于主方程的算法的计算思路为：将总的结合自由能分解为不同的能量项的贡献，并且各个能量项之间互不交叉，通过计算各能量项分量的贡献就可以得到总的结合自由能。与传统的自由能计算方法FEP和TI相比，这种方法最大的优势是可以计算配体和受体结合过程的结合自由能，而不仅仅是算出两个配体结合的相对结合自由能，并且该方法的计算量相对较小。1998年Srinivasan等在上述方法的基础上提出了基于分子力学与连续介质模型的自由能计算方法，即MM/PBSA方法。这种方法最初只用来研究核酸分子的结构稳定性，后来发展到配体-受体、蛋白质-蛋白质以及蛋白质-核酸等生物大分子体系的结合自由能计算，并得到广泛的应用。下面将对MM/PBSA方法的操作步骤进行简要的介绍。
采用MM/PBSA对体系进行结合自由能计算，以及氨基酸残基能量分解前，需要对体系进行一定时间的分子动力学模拟。结合自由能计算主要分为5个步骤，分别为计算得到轨迹、计算复合物的分子稳定性、计算结合自由能、计算熵、残基能量分解，操作步骤如下：
第一步是在分子动力学模拟之后得到复合物体系不含水的运动轨迹complex_nowat.mdcrd的基础上，通过分析轨迹，并按照一定的频率采集构象数据。在Linux界面下创建文件夹01_GenerateSnapshots，在该文件夹下执行命令mm_pbsa.pl mm_pbsa.in > mm_pbsa.log，其中mm_pbsa.in的代码内容如下：
图2 自由能计算第一步mm_pbsa.in文件的内容
其中，各参数解释为PREFIX是保存文件的名称；COMPT，RECPT，LIGPT分别表示复合物，受体蛋白和配体小分子的拓扑文件的绝对路径（查看文件绝对路径的方法为，在该文件目录下输入pwd命令即可）；GC表示该步骤的快照是从分子模拟的轨迹生成的；AS表示在生成轨迹快照的过程中发生了氨基酸残基的突变；DC表示将自由能分解为不同项目的贡献；MM表示采用Amber中的sander程序计算体系在气相中的能量；GB表示在sander程序采用GB模型计算体系的去溶剂化自由能；PB表示采用delphi方法计算去溶剂化自由能，并采用pbsa中的NPOPT选项计算非极性溶剂化能；MS表示使用molsurf计算非极性溶剂化能的贡献；NM是指采用nmode或NAB方法计算体系的熵。以上参数属于计算中最重要的参数设定，即选择每一步的计算方法，1和0分别表示肯定和否定。
BOX中的NTOTAL表示复合物体系的总原子数（不含水分子，原子数可以在拓扑文件中查看）；NSTART和NSTOP分别表示进行结合自由能计算的MD模拟的时间段，6000-10000表示将6000 ps到10000 ps的轨迹进行自由能计算（具体时间段需要根据实际情况修改，一般来讲，以体系较为稳定的时间段为最佳）；NFREQ表示选取构象的频率，250表示每250步保存一次构象数据用于自由能计算（因此6000-10000步一共打印16个构象数据）；NUMBER_LIG_GROUPS中的LSTART和LSTOP分别表示配体小分子在复合物体系中起始和结束的原子序号；NUMBER_REC_GROUPS中的RSTART和RSTOP分别表示不含配体的受体蛋白起始和结束的原子序号（受体一般是从1开始，体系的原子数以配体作为结束）；TRAJECTORY的内容表示进行结合自由能计算参考的复合物的MD模拟轨迹文件的绝对路径（该轨迹文件必须与前面的拓扑文件相对应）；
第二步是在第一步计算得到体系的轨迹快照的基础上，计算体系分子稳定性。首先在Linux操作界面下建立02_MMPBSA_Stability文件夹，运行mm_pbsa.pl mm_pbsa.in > mm_pbsa.log，mm_pbsa.in文件内容如下：
图3 自由能计算第二步mm_pbsa.in文件的内容
其中，前面部分的内容与第一步一致，只有参数数值的变化。PROC表示采用哪一种方法来解PB方程，Amber默认值为2；REFE表示选取哪一种参考状态来计算PB方程，默认值为0；INDI为溶质的介电常数；EXDI表示周围溶剂的介电常数；SCALE表示格子的间距；LINIT表示计算线性PB方程的迭代次数；PRBRAD表示溶剂的探针半径；RADIOPT表示设置PB方程计算的半径，0表示从前面的拓扑文件中选择半径；SURFTEN和SURFOFF表示计算非极性自由能贡献的值；DIELC表示静电相互作用的介电常数；IGB表示选择Tsui（1）或Onufriev（2，5）的GB方法；GBSA表示选择LCPO（1）或ICOSA（2）方法计算溶剂可接近表面积；SALTCON表示溶液中可移动的抗平衡离子的浓度（M）；EXTDIEL、INTDIEL与前面的EXTD和INTD相对应；PROBE表示用于计算SAS的探针半径；DELPHI表示用于计算PB方程的delphi程序的路径。以上内容中的距离单位均为Å。
