精确计算绝对结合自由能

个人认为这是一篇非常经典的结合自由能计算的实例,并且还都有原始参数数据,可以让我们更好的了解具体设置以及注意事项。因此我会花较多的精力去了解该文章,翻译并讲解体系参数。

原文翻译

原始文献DOI: 10.1039/c5sc02678d

精确计算药物分子的绝对结合自由能

摘要: 数十年来,精确预测结合亲和度一直是计算化学的核心目标,但至今仍难以实现。尽管取得了良好的进展,但对于类药物分子来说,在药物发现中使用所需要的准确性并没有达到。在这里,我们基于热力学循环对一组不同的抑制剂结合到含溴结构域蛋白4(BRD4)进行绝对自由能计算,并证明平均绝对误差值$0.6\;kcal\;mol^{-1}$能够实现。我们也展示了类似的精确度($1.0\;kcal\;mol^{-1}$)可以在伪前沿方法中实现。溴结构域是一种表观遗传标记读取器,可以识别乙酰化基序并调节基因转录,目前正作为癌症和炎症的治疗靶点进行研究。这种前所未有的准确性为药物类化合物结合自由能的预测提供了令人兴奋的前景。

介绍:计算药物设计的“圣杯”之一是准确预测药物与其靶蛋白的亲和力。尽管药物活性分子的发展是一个多因子的优化问题,其他方面的考虑,如生物利用度和毒性发挥着重要的作用,高亲和力的化合物其预期生物目标是一个必要的要求实现一个强有力的,有选择性的,最终有效的药物。不幸的是,即使在结构信息可用的情况下,溶剂效应、蛋白质或配体构象的变化以及熵焓补偿使配体与大分子缔合过程的合理化成为一项非常复杂的任务。然而,由于理论和计算的重要进展,特别是在过去的十年中,使用基于物理的计算机模拟来预测结合亲和度,通过通过自然地考虑由于溶剂的离散性和结合时熵的变化而产生的复杂效应,是有希望实现可靠的结合能估计。
 在显式溶剂中,基于全原子分子动力学(MD)模拟的炼金术自由能计算和导向方法是在最高的理论严密性水平上运行的典型方法,并且也可用于当前典型水平的计算能力。“炼金术”方法通常又指自由能扰动(FEP),是基于非物理热力学循环,其中结合自由能被计算为配体从不同的环境(如结合态和非结合态)“插入”或“移除”多个步骤的和。通过施加一种力将配体从蛋白质中拉出,导向或拉拽法的方法遵循的是一条物理路径。这通常是通过使用Jarzynski关系的非平衡模拟来实现的,或者通过简谐地限制配体在离结合腔不同距离处的运动,然后计算平均势。其他流行的方法包括端点方法,包括显式溶剂模拟的隐式溶剂后处理,例如使用Poisson–Boltzmann或广义玻恩表面积(MM/PBSA和MM/GBSA)方法。另一种很有前途的方法是metadynamics,它具有漏斗形的抑制电位,其中添加了偏置能量,以采样多个结合情况。
 用炼金术方法计算了几种蛋白质-配体体系的绝对结合自由能。T4溶菌酶工程结合口袋是目前研究最多的大分子体系之一。Mobley et al.研究了13种单环类碎片配体与L99A疏水性T4溶菌酶空腔突变体的结合,得到了与等温滴定量热法(ITC)相比的均方根(RMS)误差大约1.9 kcal/mol。Boyce et al.研究了类似碎片状配体与T4溶菌酶微极性模型腔的结合(L99A/M102Q),与ITC得到的五种化合物具有可测亲和力相比,获得RMS误差约为1.8 kcal/mol,。另一个流行的测试系统是FK506-结合蛋白质(EKBP12)。Holt et al.用FKBP12对一系列配体进行了初步的实验研究,和一样是药物,有多个环和几个可旋转的键,尽管它们具有非常相似的化学成分。对于九种抑制剂的亲和预测,Shirts最先报道的RMS误差大约为2.0 kcal/mol,紧接着是Wang et al.得到2.0-2.5 kcal/mol误差。该体系以竞争抑制FKPB12活性的实验自由能作为参考。Fujitani和他的同事得到的八种FKBP-12抑制剂的线性拟合RMS差异仅为0.4 kcal/mol,然而相对于实验有个大的偏离(-3.2 kcal/mol)。其他的计算也被报道过,尽管是在更少的配体上,这使得建立实际误差变得更加困难。
 兴趣驱动支持化学工具通过模拟方法的发展,我们寻求评估是否基于标准实现的炼金术转换的绝对自由能计算现在正在达到的地步可以应用于多样化,药物类有机分子和药物相关目标。为了达到我们的目标,我们比较了11种不同的小分子抑制剂与溴结构域(BRDs)结合的预测结合自由能与主要是等温滴定量热法(ITC)的实验测量结果。

BRDs是表观遗传标记读取器,能够特定的识别$\varepsilon$-N-lysine乙酰化图案(图1),已在46种人类细胞核和细胞质蛋白中发现。乙酰化作用通常发现参与染色质重塑的大分子复合物,DNA修复和细胞周期控制,尤其是在组蛋白。组蛋白乙酰化被认为会导致转录激活,乙酰化水平的改变与癌症的异常转录和呼吸困难有关。因此新型BRD抑制剂正广泛应用于医学和基础生物学研究,事实上,各种BRD抑制剂目前正处于一期和二期临床试验,用于治疗NUT中线癌、急性白血病、进展性淋巴瘤和动脉硬化。
 在此,我们进行回顾性分析,以评估在最佳情况下计算的准确性和精密度方面的性能。随后,我们进行了一项伪前瞻性研究,在不使用任何蛋白质配体复合物的结构信息的情况下,用传统的对接方法重复练习,给出初始姿势。两项研究均与实验数据吻合良好。我们讨论的结果,就如何使用这样的计算可以帮助药物的发现和开发过程。

方法:系统设置 除了11个配体外,配合物的初始构象来自于完整晶体结构(3U5J, 3U5L, 4OGI, 4OGJ, 3MXF, 4MR3, 4MR4, 3SVG, 4J0R, 4HBV),是基于结构3SVG建模,来自配体对接到apo蛋白(2OSS))的结果。晶体中缺失的原子是使用WHAT-IF网页接口建模的,所有非配体的有机分子都从体系中除去,而所有的晶体水都被保留下来。然后通过Maestro (v9.5, Schr¨odinger)添加氢,配体参数来自通用的AMBERl力场(GAFF v1.5),AM1-BCC电荷使用AmberTools12提供的FESetup tool v1.1pre1(http://ccpforge.cse.rl.ac.uk/gf/project/ccpbiosim) 得到。GROMACS拓扑和坐标来自AMBER,使用acpype(v.2013-11-28 Rev: 399)。对于蛋白和TIP3P水模型分子我们使用Amber99SB-ILDN力场。复合物溶解在十二面体盒子中,除了在立方盒中溶解的配体1和4。使用周期性边界条件,溶质分子和盒子边缘的最小距离为12埃米。加入钠离子和氯离子使体系成电中性。
自由能计算 利用图2所示的非物理热力学循环,分别从晶体配体姿态和对接姿态开始进行绝对结合自由能的计算。

所有的模拟都在GROMACS 4.6.5中进行。通过使用线性炼金术途径,对范德华相互作用使用$\Delta\lambda = 0.05$,库伦相互作用使用$\Delta\lambda = 0.1$,使得配体范德华相互作用被解耦合,电荷湮灭。为了取代配体约束,我们使用了12个非均匀分布的$\lambda$值(0.0, 0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.5, 0.75, 1.0)。因此,总共有42个窗口用于复合物模拟,31个窗口用于配体模拟。对于每个窗口,使用最速下降算法执行10000能量最小化步骤。然后在施加力常数为1000 kJ/mol/nm的溶质重原子的简谐位置约束的正则系综中,对系统进行了0.5 ns的数值模拟。采用Langevin动力学,298.15 K作为参考温度进行温度耦合。在等温-等压系综中,采用Berendsen弱耦合算法进行了1 ns位置约束运行。采用Hamiltonian-exchange Langevin动力学方法,以2 fs时间步长,在NPT系综下进行了10 ns无约束成品运行用于数据采集。每1000次尝试在任意状态对之间进行300万次交换,按照Chodera和Shirts等人提出的吉布斯抽样方案。这使得相邻状态的接收率从0.1到0.3不等(均值和标准差为$0.2\pm0.1$),但是跳跃到其他状态的概率从0.2到0.9(均值和标准差为$0.7\pm0.2$)。因此,本研究的总的无限制模拟时间为43 ms。结合配体相对于蛋白的相对位置和方向通过一个距离,两个角度和三个二面体简谐势力常数为10 kcal/mol/A^(-2)来限制。这组约束对自由能的贡献可以按照Boresch et al.描述的对于溶液中的非相互作用配体($\Delta G^{solv}_{restr}$)解析的方法计算出来,而对于与蛋白质相互作用的配体复合物则需要进行数值评估。用于评估这种贡献的方程还包括对结合自由能的标准态依赖性的校正。一种”软核”势用于范德瓦尔斯相互作用变换。在所有的模拟中,粒子网格Ewald (PME)算法用于静电相互作用,实体空间截断为12埃米。一个6阶的样条,相对容差为$10^{-6}$,傅里叶间距1.0埃米。9-10埃米的switch函数用于范德华相互作用。P-LINCS约束算法仅仅用于H-bonds。按照Shirts et al,使用了能量和压力的GROMACS远程色散校正,以及附加的远程色散校正(EXP-LR)。对于带点化合物(1和4),采用Rocklin et al.提出了静电场尺寸效应的半解析校正方案。使用APBS 1.3计算了11帧的剩余综合电位(RIP),平均RIP值仅用于校正;由于配位体和络合物的构象不同,其标准偏差始终在以下0.05 kcal/mol(每1ns,从1ns到10ns)。
数据分析 采用多重Bennet acceptance ratio(MBAR)对结果进行了分析,由python包提供(https://simtk.org/home/pymbar) 。每个窗口的前 1 ns被丢弃。在自由能估计之前,通过计算势能的自相关时间和统计效率,对每个状态的数据进行子采样,只包含不相关的数据点。采用随机重采样的方法构建了200个boot-strap set,并替换了不相关的数据,其中第一个set是原始样本。用MBAR法计算了所有200组的自由能,最终的估计值是所有这些自由能估计值的平均值。对于所有的bootstrap采样,最终的误差估计是其采样的标准偏差。对于从晶体配体姿态开始的计算,整个计算重复了三次,以评估结果的收敛性。利用600个bootstrap样本的均值和样本标准差,确定配体和复合物模拟的自由能估计值和标准差。因此,最终的不确定度是MBAR自由能估计的统计不确定性和由有限采样引起的误差的两者表示。配体11是个例外,因为它是基于配体9和PDB结构3SVG建模的,存在两种同样合理的结合模式,这儿三氟甲苯部分翻转180度。因此,每一种结合模型进行两次计算重复,共得到四次结合自由能计算。根据Mobley et al.描述。多种结合模式的结果可以合并为一个单一的结合自由能值。配体与水溶液的解耦以及与溶剂化络合物的解耦之间的对每个配体的最终结合自有能是有差异的。因此,结合自由能的最终误差是配体和复合物计算的不确定性的平方根。
对接 使用MOE v2013.08对apo蛋白(BRD4(1))的一个晶体结构(PDB 2OSS)进行刚性对接。存在于所有BET溴域结合口袋的五个高度保守的结晶水被保留,而其他所有的水和有机分子都被去除。配体的二维结构由Marvin Sketch (v6.1.0, ChemAxon)绘制,为了生成三维构象进行了随机构象搜索。每个配体的构象最多限制为100个,去除重复构象(RMSD < 0.25埃米)。对接方案采用药包放置方法和$\Delta G$打分函数。然后将每个结合姿势最小化并使用GBVI/WSA$\Delta G$打分函数进行重打分。药效团查询是基于PDB结构3UVW中发现的乙酰赖氨酸的性质建立的,由一个模拟乙酰氧的氢键受体位点和一个与甲基半基位置对应的非极性位点组成。蛋白质使用Amber ff99SB进行参数化。利用二维扩展H¨uckel理论得到配体键合参数,根据MOE中的AMBER10:EHT力场选项,VdW参数来自GAFF,电荷来自键电荷增量。重复的姿态会根据它们的氢键和疏水模式自动移除,GBVI/WSA$\Delta G$打分函数预测的具有正结合能的姿态同样被移除,因为它们通常与蛋白质原子发生明显的冲突。其余的姿态进一步由RMSD以3埃米截断进行聚类,以获得可能的结合姿态的大体形貌。同时考虑到我们对在模拟时间尺度内可以相互转换的类似结合方向上的自由能计算不感兴趣。因此,只有在每个集群内的最佳得分姿态将用于自由能计算。这一过程旨在减少计算量,同时最大限度地保留最接近晶体的姿态。

结果 基于晶体结构的绝对自由能计算是精确的 在这项研究中,我们使用非物理热力学循环进行了彻底的结合自由能计算(图2),从晶体结构BRD4(1)开始,与11中抑制剂形成了相同的结合位点(图1,ESI图1)。我们最先在有利情况下对评估预测性能感兴趣,也就是说,当结合构象从实验中已知。因此,这项研究的结果提供了一幅可以预期的最佳准确度的图像。另外,我们对评估计算的精度很感兴趣,考虑到在测试集中存在大而易溶的分子(图3)。在处理这类药物分子时,结果中存在很大的不确定性,这确实会妨碍对结果的准确性进行有意义的评估。为此,除了bootstrap分析去评价自由能估计量的统计不确定性,我们决定重复计算三次,以便得到不确定度的近似值(如方法部分所述,对配体11进行了四次重复)。事实上,虽然bootstrap提供了比单独的MBAR误差估计更现实的不确定性估计,但它仍然低估了样本标准差。每次计算结果为73,10 ns长,这部分研究全原子分子动力学运行总模拟时间约为25微秒。
 考虑的抑制剂主要包括具有多种理化性质的类药物分子:原子数从22到77;分子质量重241到525 Da;可旋转键数从0到11;计算的logP从-0.4到5.3(ESI表1)。亲和力的范围包括微摩尔的结合,例如配体10(~23 uM)和11(~80 uM);低到纳米级的结合,例如配体1(~40 nM)和2(~50 nM)。许多不同的化学基团被表示出来,而集合的不同提供给我们更确信所得到的结果不会因考虑到有限的化学空间而产生过大的偏差。
 表1总结了这一回顾性研究的结果(参见ESI表2对能量贡献的分解)。

大多数计算结果与实验测定值非常吻合。11项预测中有7项误差在0.5 kcal/mol以下,所有预测误差都低于2.0 kcal/mol。这导致了一个平均绝对误差(MAE)$0.6\pm0.1 kcal/mol$,均方根误差(RMS)$0.8\pm0.2 kcal/mol$。计算得到的自由能与实验值有较强的相关性如表1和图5所示,Pearson指数r为$0.84\pm0.05$,并对配体亲和度进行有效排序(Spearman $p=0.82\pm0.06$)。计算的精确度也令人鼓舞,因为只有三种情况的不确定性高于上述0.5 kcal/mol水平,且所有情况都低于1.0 kcal/mol(见ESI图2)。正如所料,最大的不确定性发生在考虑最大配体的时候。
在缺乏复杂结构信息的情况下,也可以实现准确的预测
 为了评估此类计算在未来环境中的实用性,我们进行了基于对接结果而不是基于晶体结构的对接和自由能计算练习。这部分研究的主要目的是评估在回顾性研究中观察到的准确性是否可以通过更前瞻性的基于对接的方法来实现。三维配体结构是通过构象搜索在二维空间中绘制配体而产生的。然后创建了一个基于乙酰赖氨酸的药效团,用于帮助11种抑制剂与BRD4(1)(pdb-ID 2OSS)的apo结构对接。无约束分子动力学仿真研究的是紧密结合的构象,因此得到的对接姿态采用均方根偏差(RMSD)聚类,以避免在仿真时间尺度内选择相互转换的构象。所获得的对接姿态总数为72个,根据聚类过程进行分析减少,根据评分函数的预测去除具有正结合自由能的姿态。图4显示了这25个对接姿态(crystal得分和RMSD总结在ESI表3中)。

结论 我们在这里展示了,对于一个小的,相当刚性的系统,比如溴结合域,基于分子动力学的自由能计算能够实现RMS误差当启始结构来自对接时不超过1.4 kcal/mol,当使用晶体结构和更昂贵的方案时下降到0.8 kcal/mol。目前的结果证实了绝对自由能计算在药物发现应用中的潜力。据我们所知,这是第一个关于绝对结合自由能的研究,该研究考虑了一系列不同的类药物分子和目前研究其治疗潜力的生物学相关靶点。值得注意的是,最近有报道称,在相对结合自由能方面,对大量分子也有类似的精确度。绝对自由能计算的可靠性保证了它们至少在相当严格的药物靶标(如溴结构域)的药物发现运动中得到应用。

参考文献 略,见原始文献。

原始数据解析

数据源结构如下,主要包括4个部分的文件:

  • crystal—包含计算所需的所有拓扑和坐标文件,即以复合物的结构为起始结构。
  • docking—包含计算所需的所有拓扑和坐标文件,这是基于对接产生的起始结构。
  • mdp—包含用于所有平衡和生产运行的.mdp文件,无论是溶液中的配体还是复合物。
  • xvg—包含压缩大小的dhdl。通过模拟得到了xvg文件。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
├─crystal
│ ├─ligand1
│ │ ├─complex
│ │ │ └─example_xvg_files
│ │ └─ligand
│ │ └─example_xvg_files
│ ├─Ligand10
│ │ ├─complex
│ │ └─ligand
│ ├─ligand11
│ │ ├─model1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─Model2
│ │ ├─complex
│ │ └─ligand
│ ├─Ligand2
│ │ ├─complex
│ │ └─ligand
│ ├─ligand3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand4
│ │ ├─complex
│ │ └─ligand
│ ├─Ligand5
│ │ ├─complex
│ │ └─ligand
│ ├─ligand6
│ │ ├─complex
│ │ └─ligand
│ ├─Ligand7
│ │ ├─complex
│ │ └─ligand
│ ├─Ligand8
│ │ ├─complex
│ │ └─ligand
│ └─Ligand9
│ ├─complex
│ └─ligand
├─docking
│ ├─ligand1
│ │ ├─pose1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─pose2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─Pose3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand10
│ │ └─pose1
│ │ ├─complex
│ │ └─ligand
│ ├─Ligand11
│ │ ├─Pose1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─Pose2
│ │ ├─complex
│ │ └─ligand
│ ├─Ligand2
│ │ ├─Pose1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─Pose2
│ │ ├─complex
│ │ └─ligand
│ ├─ligand3
│ │ ├─Pose1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─pose2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─pose3
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─Pose4
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─pose5
│ │ ├─complex
│ │ └─ligand
│ ├─Ligand4
│ │ └─Pose1
│ │ ├─complex
│ │ └─ligand
│ ├─ligand5
│ │ ├─Pose1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─pose2
│ │ ├─complex
│ │ └─ligand
│ ├─ligand6
│ │ ├─Pose1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─pose2
│ │ ├─complex
│ │ └─ligand
│ ├─ligand7
│ │ ├─pose1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─Pose2
│ │ ├─complex
│ │ └─ligand
│ ├─ligand8
│ │ ├─pose1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─Pose2
│ │ ├─complex
│ │ └─ligand
│ └─ligand9
│ ├─Pose1
│ │ ├─complex
│ │ └─ligand
│ ├─pose2
│ │ ├─complex
│ │ └─ligand
│ └─pose3
│ ├─complex
│ └─ligand
├─mdp
│ ├─complex
│ │ ├─ENMIN
│ │ ├─NPT
│ │ ├─NVT
│ │ └─PROD
│ └─ligand
│ ├─ENMIN
│ ├─NPT
│ ├─NVT
│ └─PROD
└─xvg
├─crystal
│ ├─ligand10_quinaz
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand11_isoxcf3
│ │ ├─model1
│ │ │ ├─run1
│ │ │ │ ├─complex
│ │ │ │ └─ligand
│ │ │ └─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─Model2
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run2
│ │ ├─complex
│ │ └─ligand
│ ├─ligand1_bi2356
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand2_jq1
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ ├─Ligand3_RVXOH
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand4_tg101348
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand5_isox2
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand6_bzt7
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand7_rvx208
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ ├─ligand8_alprazolam
│ │ ├─run1
│ │ │ ├─complex
│ │ │ └─ligand
│ │ ├─run2
│ │ │ ├─complex
│ │ │ └─ligand
│ │ └─run3
│ │ ├─complex
│ │ └─ligand
│ └─ligand9_isox1
│ ├─run1
│ │ ├─complex
│ │ └─ligand
│ ├─run2
│ │ ├─complex
│ │ └─ligand
│ └─run3
│ ├─complex
│ └─ligand
└─docking
├─Ligand10_Quinaz
│ └─Pose1
│ ├─complex
│ └─ligand
├─Ligand11_IsoxCF3
│ ├─Pose1
│ │ ├─complex
│ │ └─ligand
│ └─pose2
│ ├─complex
│ └─ligand
├─Ligand1_BI2356
│ ├─Pose1
│ │ ├─complex
│ │ └─ligand
│ ├─pose2
│ │ ├─complex
│ │ └─ligand
│ └─pose3
│ ├─complex
│ └─ligand
├─ligand2_jq1
│ ├─Pose1
│ │ ├─complex
│ │ └─ligand
│ └─Pose2
│ ├─complex
│ └─ligand
├─Ligand3_RVXOH
│ ├─Pose1
│ │ ├─complex
│ │ └─ligand
│ ├─Pose2
│ │ ├─complex
│ │ └─ligand
│ ├─Pose3
│ │ ├─complex
│ │ └─ligand
│ ├─pose4
│ │ ├─complex
│ │ └─ligand
│ └─Pose5
│ ├─complex
│ └─ligand
├─Ligand4_TG101348
│ └─Pose1
│ ├─complex
│ └─ligand
├─ligand5_isox2
│ ├─pose1
│ │ ├─complex
│ │ └─ligand
│ └─Pose2
│ ├─complex
│ └─ligand
├─Ligand6_BzT7
│ ├─pose1
│ │ ├─complex
│ │ └─ligand
│ └─Pose2
│ ├─complex
│ └─ligand
├─Ligand7_RVX208
│ ├─pose1
│ │ ├─complex
│ │ └─ligand
│ └─Pose2
│ ├─complex
│ └─ligand
├─Ligand8_Alprazolam
│ ├─pose1
│ │ ├─complex
│ │ └─ligand
│ └─Pose2
│ ├─complex
│ └─ligand
└─Ligand9_Isox1
├─Pose1
│ ├─complex
│ └─ligand
├─Pose2
│ ├─complex
│ └─ligand
└─pose3
├─complex
└─ligand

前两项就不说了,就是初始结构。下面主要讲解后面两项在GROMACS模拟中用到的文件。

mdp参数文件

数据结构如下:

1
2
3
4
5
6
7
8
9
10
├─complex
│ ├─ENMIN
│ ├─NVT
│ ├─NPT
│ └─PROD
└─ligand
├─ENMIN
├─NVT
├─NPT
└─PROD

不难看出是分别对复合物和配体进行能量极小化、NVT预平衡、NPT预平衡以及成品模拟。下面我们一个个来看看参数设置:

复合物

  • 能量极小化mdp,相关参数说明我直接附在后面了。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    ;====================================================
    ; Energy minimization
    ;====================================================

    ; RUN CONTROL & MINIMIZATION
    ;----------------------------------------------------
    define = -DFLEXIBLE ; 使用柔性水模型
    integrator = steep ; 采用最速下降法
    nsteps = 10000 ; 最大步数
    emtol = 10 ; 容忍值需要达到10 kJ/mol
    emstep = 0.01 ; 步长
    nstcomm = 100 ; 移除质心运行的频率(步)

    ; OUTPUT CONTROL
    ;----------------------------------------------------
    nstxout = 200 ; 每200步保存坐标到.trr文件中
    nstvout = 0 ; 不需要保存速度
    nstfout = 0 ; 不需要保存力

    nstxtcout = 500 ; 每500步输出压缩轨迹.xtc
    xtc-precision = 1000 ; 压缩轨迹精度
    nstlog = 500 ; 每500更新日志文件
    nstenergy = 500 ; 每500步输出能量文件
    nstcalcenergy = 100 ; 计算能量的频率

    ; NEIGHBOR SEARCHING
    ;----------------------------------------------------
    cutoff-scheme = group ; 'group'在gromacs4.6默认,groamcs5以上默认为 'Verlet'
    ns-type = grid ; 邻区搜索算法
    nstlist = 1 ; 邻区列表更新频率
    rlist = 1.2 ; 邻区列表截断距离(nm)

    ; BONDS
    ;----------------------------------------------------
    constraints = none ; 仅拓扑中指定的约束

    ; ELECTROSTATICS
    ;----------------------------------------------------
    coulombtype = PME ; 静电计算方法
    rcoulomb = 1.2 ; 静电截断半径
    pme-order = 6 ; 内插阶数
    fourierspacing = 0.10 ; PME/PPPM方法FFT格点最大间距(nm)
    ewald-rtol = 1e-6 ; 静电能量容差

    ; VDW
    ;----------------------------------------------------
    vdw-type = Switch ; 范德华计算方法
    rvdw-switch = 0.9 ; 切换距离
    rvdw = 1.0 ; 范德华截断半径
    DispCorr = EnerPres ; 长程色散校正, 能量和压力

    ; TEMPERATURE & PRESSURE COUPL
    ;----------------------------------------------------
    Tcoupl = no ; 耦合方法, 无
    Pcoupl = no ; 压力耦合, 无
    gen_vel = no ; 不产生速度

    ; FREE ENERGY
    ;----------------------------------------------------
    free-energy = yes ; 进行自由能计算
    sc-alpha = 0.5 ; 软核势参数
    sc-power = 1 ; 次数1
    sc-sigma = 0.3 ; 用于C6/C12为零
    sc-coul = no ; 对静电不使用软核
    init-lambda-state = 0 ; 从零开始的状态编号
    restraint-lambdas = 0.0 0.01 0.025 0.05 0.075 0.1 0.15 0.2 0.35 0.5 0.75 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 ; 限制
    coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 ; 静电作用
    vdw-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 ; vdw作用
    nstdhdl = 100 ; 输出频率, 为nstcalcenergy的倍数
    dhdl-print-energy = yes ; 包含总能量或势能
    calc-lambda-neighbors = -1 ; 相邻的数目. mbar取-1
  • NVT预平衡mdp:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90

    ;====================================================
    ; NVT equilibration
    ;====================================================

    ; RUN CONTROL
    ;----------------------------------------------------
    define = -DPOSRES ; 位置限制
    integrator = sd ; 随机动力学
    nsteps = 250000 ; 总步数,共250000*0.002 = 500 ps
    dt = 0.002 ; 步长2 fs
    comm-mode = Linear ; 移除质心运动的方式, Linear: 平动
    nstcomm = 100 ; 移除质心运动频率

    ; OUTPUT CONTROL
    ;----------------------------------------------------
    nstxout = 50000 ; 每50000步保存轨迹到.trr
    nstvout = 0 ; 不保存速度
    nstfout = 0 ; 不保存力

    nstxtcout = 1000 ; 压缩轨迹输出频率
    xtc-precision = 1000 ; 压缩轨迹输出精度
    nstlog = 1000 ; 更新log日志频
    nstenergy = 1000 ; 保存能量频率
    nstcalcenergy = 100 ; 计算能量的频率

    ; BONDS
    ;----------------------------------------------------
    constraint_algorithm = lincs ; 约束算法, lincs: 不支持键角
    constraints = h-bonds ; 仅氢键约束
    lincs_iter = 1 ; lincs最终步的迭代数, 1: 常规模拟
    lincs_order = 6 ; 约束耦合矩阵展开的最高阶数
    lincs-warnangle = 30 ; 如果某一步中键的旋转角度超过此值, 显示lincs警告
    continuation = no ; 对初始构象应用约束和重置shell

    ; NEIGHBOR SEARCHING
    ;----------------------------------------------------
    cutoff-scheme = group ; 'group'在gromacs4.6默认,groamcs5以上默认为 'Verlet'
    ns-type = grid ; 邻区搜索算法
    nstlist = 10 ; 邻区列表更新频率
    rlist = 1.2 ; 邻区列表截断距离(nm)
    pbc = xyz ; 周期性边界条件, xyz;

    ; ELECTROSTATICS & EWALD
    ;----------------------------------------------------
    coulombtype = PME ; 静电计算方法
    rcoulomb = 1.2 ; 静电截断半径
    ewald_geometry = 3d ; Ewald几何结构
    pme-order = 6 ; 内插阶数
    fourierspacing = 0.10 ; PME/PPPM方法FFT格点最大间距(nm)
    ewald-rtol = 1e-6 ; 静电能量容差
    optimize-fft = no ; 不要在启动时计算网格的最优FFT计划

    ; VAN DER WAALS
    ;----------------------------------------------------
    vdwtype = Switch ; 范德华计算方法
    rvdw = 1.0 ; 范德华截断半径
    rvdw-switch = 0.9 ; 切换距离
    DispCorr = EnerPres ; 长程色散校正, 能量和压力

    ; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
    ;----------------------------------------------------
    tc_grps = Protein non-Protein ; 温度耦合组, 蛋白和非蛋白
    tau_t = 1.0 1.0 ; 时间常数(ps)
    ref_t = 298.15 298.15 ; 参考温度(K)

    ; PRESSURE COUPLING
    ;----------------------------------------------------
    pcoupl = no ; 无压力耦合

    ; VELOCITY GENERATION
    ;----------------------------------------------------
    gen_vel = yes ; 产生速度 (如果gen_vel是'yes', continuation应该设置为'no')
    gen_seed = -1 ; 随机数种子; -1: 自动确定
    gen_temp = 298.15 ; 随机速度对应的温度

    ; FREE ENERGY
    ;----------------------------------------------------
    free-energy = yes ; 进行自由能计算
    sc-alpha = 0.5 ; 软核势参数
    sc-power = 1 ; 次数1
    sc-sigma = 0.3 ; 用于C6/C12为零
    sc-coul = no ; 对静电不使用软核
    init-lambda-state = 0 ; 从零开始的状态编号
    restraint-lambdas = 0.0 0.01 0.025 0.05 0.075 0.1 0.15 0.2 0.35 0.5 0.75 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 ;限制
    coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 ; 静电作用
    vdw-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 ; vdw作用
    nstdhdl = 100 ; 输出频率, 为nstcalcenergy的倍数
    dhdl-print-energy = yes ; 包含总能量或势能
    calc-lambda-neighbors = -1 ; 相邻的数目. mbar取-1
  • NPT预平衡mdp

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94

    ;====================================================
    ; NPT equilibration
    ;====================================================

    ; RUN CONTROL
    ;----------------------------------------------------
    define = -DPOSRES ; 位置限制
    integrator = sd ; stochastic leap-frog integrator
    nsteps = 500000 ; 2 * 500,000 fs = 1 ns
    dt = 0.002 ; 2 fs
    comm-mode = Linear ; remove center of mass translation
    nstcomm = 100 ; frequency for center of mass motion removal

    ; OUTPUT CONTROL
    ;----------------------------------------------------
    nstxout = 50000 ; save coordinates to .trr every 100 ps
    nstvout = 0 ; don't save velocities to .trr
    nstfout = 0 ; don't save forces to .trr

    nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps
    xtc-precision = 1000 ; precision with which to write to the compressed trajectory file
    nstlog = 1000 ; update log file every 2 ps
    nstenergy = 1000 ; save energies every 2 ps
    nstcalcenergy = 100 ; calculate energies every 100 steps

    ; BONDS
    ;----------------------------------------------------
    constraint_algorithm = lincs ; holonomic constraints
    constraints = h-bonds ; hydrogens only are constrained
    lincs_iter = 1 ; accuracy of LINCS (1 is default)
    lincs_order = 6 ; also related to accuracy (4 is default)
    lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
    continuation = no ; 对初始构象应用约束和重置shell ?

    ; NEIGHBOR SEARCHING
    ;----------------------------------------------------
    cutoff-scheme = group ; 'group' is default in 4.6 - from GMX5 on 'Verlet' is default
    ns-type = grid ; search neighboring grid cells
    nstlist = 10 ; 20 fs (default is 10)
    rlist = 1.2 ; short-range neighborlist cutoff (in nm)
    pbc = xyz ; 3D PBC

    ; ELECTROSTATICS & EWALD
    ;----------------------------------------------------
    coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
    rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
    ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
    pme-order = 6 ; interpolation order for PME (default is 4)
    fourierspacing = 0.10 ; grid spacing for FFT
    ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb
    optimize-fft = no ; don't calculate the optimal FFT plan for the grid at startup

    ; VAN DER WAALS
    ;----------------------------------------------------
    vdwtype = Switch
    rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
    rvdw-switch = 0.9
    DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure

    ; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
    ;----------------------------------------------------
    tc_grps = Protein non-Protein
    tau_t = 1.0 1.0
    ref_t = 298.15 298.15

    ; PRESSURE COUPLING
    ;----------------------------------------------------
    pcoupl = Berendsen ; 压力耦合方法采用Berendsen
    pcoupltype = isotropic ; 耦合类型, isotropic: 各向同性;
    tau_p = 0.5 ; 时间常数(ps)
    ref_p = 1.0 ; 参考压力(bar)
    compressibility = 4.5e-05 ; 压缩率(1/bar)

    ; VELOCITY GENERATION
    ;----------------------------------------------------
    gen_vel = yes ; Velocity generation is on (if gen_vel is 'yes', continuation should be 'no') ?
    gen_seed = -1 ; Use random seed
    gen_temp = 298.15

    ; FREE ENERGY
    ;----------------------------------------------------
    free-energy = yes
    sc-alpha = 0.5
    sc-power = 1
    sc-sigma = 0.3
    sc-coul = no
    init-lambda-state = 0
    restraint-lambdas = 0.0 0.01 0.025 0.05 0.075 0.1 0.15 0.2 0.35 0.5 0.75 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
    coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
    vdw-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
    nstdhdl = 100
    dhdl-print-energy = yes
    calc-lambda-neighbors = -1
  • 成品模拟mdp

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93

    ;====================================================
    ; Production MD
    ;====================================================

    ; RUN CONTROL
    ;----------------------------------------------------
    integrator = sd ; stochastic leap-frog integrator
    nsteps = 5000000 ; 2 * 5,000,000 fs = 10 ns
    dt = 0.002 ; 2 fs
    comm-mode = Linear ; remove center of mass translation
    nstcomm = 100 ; frequency for center of mass motion removal

    ; OUTPUT CONTROL
    ;----------------------------------------------------
    nstxout = 50000 ; save coordinates to .trr every 100 ps
    nstvout = 0 ; don't save velocities to .trr
    nstfout = 0 ; don't save forces to .trr

    nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps
    xtc-precision = 1000 ; precision with which to write to the compressed trajectory file
    nstlog = 1000 ; update log file every 2 ps
    nstenergy = 1000 ; save energies every 2 ps
    nstcalcenergy = 100 ; calculate energies every 100 steps

    ; BONDS
    ;----------------------------------------------------
    constraint_algorithm = lincs ; holonomic constraints
    constraints = h-bonds ; hydrogens only are constrained
    lincs_iter = 1 ; accuracy of LINCS (1 is default)
    lincs_order = 6 ; also related to accuracy (4 is default)
    lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
    continuation = no ; 对初始构象应用约束和重置shell ??

    ; NEIGHBOR SEARCHING
    ;----------------------------------------------------
    cutoff-scheme = group ; 'group' is default in 4.6 - from GMX5 on 'Verlet' is default
    ns-type = grid ; search neighboring grid cells
    nstlist = 10 ; 20 fs (default is 10)
    rlist = 1.2 ; short-range neighborlist cutoff (in nm)
    pbc = xyz ; 3D PBC

    ; ELECTROSTATICS & EWALD
    ;----------------------------------------------------
    coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
    rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
    ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
    pme-order = 6 ; interpolation order for PME (default is 4)
    fourierspacing = 0.10 ; grid spacing for FFT
    ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb
    optimize-fft = no ; don't calculate the optimal FFT plan for the grid at startup

    ; VAN DER WAALS
    ;----------------------------------------------------
    vdwtype = Switch
    rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
    rvdw-switch = 0.9
    DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure

    ; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
    ;----------------------------------------------------
    tc_grps = Protein non-Protein
    tau_t = 1.0 1.0
    ref_t = 298.15 298.15

    ; PRESSURE COUPLING
    ;----------------------------------------------------
    pcoupl = Parrinello-Rahman
    pcoupltype = isotropic ; uniform scaling of box vectors
    tau_p = 2 ; time constant (ps)
    ref_p = 1.0 ; reference pressure (bar)
    compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)

    ; VELOCITY GENERATION
    ;----------------------------------------------------
    gen_vel = yes ; Velocity generation is on (if gen_vel is 'yes', continuation should be 'no') ??
    gen_seed = -1 ; Use random seed
    gen_temp = 298.15

    ; FREE ENERGY
    ;----------------------------------------------------
    free-energy = yes
    sc-alpha = 0.5
    sc-power = 1
    sc-sigma = 0.3
    sc-coul = no
    init-lambda-state = 0
    restraint-lambdas = 0.0 0.01 0.025 0.05 0.075 0.1 0.15 0.2 0.35 0.5 0.75 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
    coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
    vdw-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
    nstdhdl = 100
    dhdl-print-energy = yes
    calc-lambda-neighbors = -1

总共包含42个状态,即init-lambda-state的值0 ~ 41,然后分别模拟即可。

配体

  • 能量极小化mdp

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71

    ;====================================================
    ; Energy minimization
    ;====================================================

    ; RUN CONTROL & MINIMIZATION
    ;----------------------------------------------------
    define = -DFLEXIBLE
    integrator = steep
    nsteps = 10000
    emtol = 10
    emstep = 0.01
    nstcomm = 100

    ; OUTPUT CONTROL
    ;----------------------------------------------------
    nstxout = 200 ; save coordinates to .trr every 200 steps
    nstvout = 0 ; don't save velocities to .trr
    nstfout = 0 ; don't save forces to .trr

    nstxtcout = 500 ; xtc compressed trajectory output every 500 steps
    xtc-precision = 1000
    nstlog = 500 ; update log file every 500 steps
    nstenergy = 500 ; save energies every 500 steps
    nstcalcenergy = 100

    ; NEIGHBOR SEARCHING
    ;----------------------------------------------------
    cutoff-scheme = group ; 'group' is default in 4.6 - from GMX5 on 'Verlet' is default
    ns-type = grid
    nstlist = 1
    rlist = 1.2

    ; BONDS
    ;----------------------------------------------------
    constraints = none

    ; ELECTROSTATICS
    ;----------------------------------------------------
    coulombtype = PME
    rcoulomb = 1.2
    pme-order = 6
    fourierspacing = 0.10
    ewald-rtol = 1e-6

    ; VDW
    ;----------------------------------------------------
    vdw-type = Switch
    rvdw-switch = 0.9
    rvdw = 1.0
    DispCorr = EnerPres

    ; TEMPERATURE & PRESSURE COUPL
    ;----------------------------------------------------
    Tcoupl = no
    Pcoupl = no
    gen_vel = no

    ; FREE ENERGY
    ;----------------------------------------------------
    free-energy = yes
    sc-alpha = 0.5
    sc-power = 1
    sc-sigma = 0.3
    sc-coul = no
    init-lambda-state = 0
    coul-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
    vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
    nstdhdl = 100
    dhdl-print-energy = yes
    calc-lambda-neighbors = -1
  • NVT预平衡mdp

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89

    ;====================================================
    ; NVT equilibration
    ;====================================================

    ; RUN CONTROL
    ;----------------------------------------------------
    define = -DPOSRES
    integrator = sd ; stochastic leap-frog integrator
    nsteps = 250000 ; 2 * 250,000 fs = 500 ps
    dt = 0.002 ; 2 fs
    comm-mode = Linear ; remove center of mass translation
    nstcomm = 100 ; frequency for center of mass motion removal

    ; OUTPUT CONTROL
    ;----------------------------------------------------
    nstxout = 50000 ; save coordinates to .trr every 100 ps
    nstvout = 0 ; don't save velocities to .trr
    nstfout = 0 ; don't save forces to .trr

    nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps
    xtc-precision = 1000 ; precision with which to write to the compressed trajectory file
    nstlog = 1000 ; update log file every 2 ps
    nstenergy = 1000 ; save energies every 2 ps
    nstcalcenergy = 100 ; calculate energies every 100 steps

    ; BONDS
    ;----------------------------------------------------
    constraint_algorithm = lincs ; holonomic constraints
    constraints = h-bonds ; hydrogens only are constrained
    lincs_iter = 1 ; accuracy of LINCS (1 is default)
    lincs_order = 6 ; also related to accuracy (4 is default)
    lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
    continuation = no ; formerly known as 'unconstrained-start' - useful for exact continuations and reruns

    ; NEIGHBOR SEARCHING
    ;----------------------------------------------------
    cutoff-scheme = group ; 'group' is default in 4.6 - from GMX5 on 'Verlet' is default
    ns-type = grid ; search neighboring grid cells
    nstlist = 10 ; 20 fs (default is 10)
    rlist = 1.2 ; short-range neighborlist cutoff (in nm)
    pbc = xyz ; 3D PBC

    ; ELECTROSTATICS & EWALD
    ;----------------------------------------------------
    coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
    rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
    ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
    pme-order = 6 ; interpolation order for PME (default is 4)
    fourierspacing = 0.10 ; grid spacing for FFT
    ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb
    optimize-fft = no ; don't calculate the optimal FFT plan for the grid at startup

    ; VAN DER WAALS
    ;----------------------------------------------------
    vdwtype = Switch
    rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
    rvdw-switch = 0.9
    DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure

    ; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
    ;----------------------------------------------------
    tc_grps = System ; 温度耦合组,系统
    tau_t = 1.0
    ref_t = 298.15

    ; PRESSURE COUPLING
    ;----------------------------------------------------
    pcoupl = no

    ; VELOCITY GENERATION
    ;----------------------------------------------------
    gen_vel = yes ; Velocity generation is on (if gen_vel is 'yes', continuation should be 'no')
    gen_seed = -1 ; Use random seed
    gen_temp = 298.15

    ; FREE ENERGY
    ;----------------------------------------------------
    free-energy = yes
    sc-alpha = 0.5
    sc-power = 1
    sc-sigma = 0.3
    sc-coul = no
    init-lambda-state = 0
    coul-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
    vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
    nstdhdl = 100
    dhdl-print-energy = yes
    calc-lambda-neighbors = -1
  • NPT预平衡mdp

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93

    ;====================================================
    ; NPT equilibration
    ;====================================================

    ; RUN CONTROL
    ;----------------------------------------------------
    define = -DPOSRES
    integrator = sd ; stochastic leap-frog integrator
    nsteps = 500000 ; 2 * 500,000 fs = 1 ns
    dt = 0.002 ; 2 fs
    comm-mode = Linear ; remove center of mass translation
    nstcomm = 100 ; frequency for center of mass motion removal

    ; OUTPUT CONTROL
    ;----------------------------------------------------
    nstxout = 50000 ; save coordinates to .trr every 100 ps
    nstvout = 0 ; don't save velocities to .trr
    nstfout = 0 ; don't save forces to .trr

    nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps
    xtc-precision = 1000 ; precision with which to write to the compressed trajectory file
    nstlog = 1000 ; update log file every 2 ps
    nstenergy = 1000 ; save energies every 2 ps
    nstcalcenergy = 100 ; calculate energies every 100 steps

    ; BONDS
    ;----------------------------------------------------
    constraint_algorithm = lincs ; holonomic constraints
    constraints = h-bonds ; hydrogens only are constrained
    lincs_iter = 1 ; accuracy of LINCS (1 is default)
    lincs_order = 6 ; also related to accuracy (4 is default)
    lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
    continuation = no ; formerly known as 'unconstrained-start' - useful for exact continuations and reruns

    ; NEIGHBOR SEARCHING
    ;----------------------------------------------------
    cutoff-scheme = group ; 'group' is default in 4.6 - from GMX5 on 'Verlet' is default
    ns-type = grid ; search neighboring grid cells
    nstlist = 10 ; 20 fs (default is 10)
    rlist = 1.2 ; short-range neighborlist cutoff (in nm)
    pbc = xyz ; 3D PBC

    ; ELECTROSTATICS & EWALD
    ;----------------------------------------------------
    coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
    rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
    ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
    pme-order = 6 ; interpolation order for PME (default is 4)
    fourierspacing = 0.10 ; grid spacing for FFT
    ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb
    optimize-fft = no ; don't calculate the optimal FFT plan for the grid at startup

    ; VAN DER WAALS
    ;----------------------------------------------------
    vdwtype = Switch
    rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
    rvdw-switch = 0.9
    DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure

    ; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
    ;----------------------------------------------------
    tc_grps = System ; 温度耦合组,系统
    tau_t = 1.0
    ref_t = 298.15

    ; PRESSURE COUPLING
    ;----------------------------------------------------
    pcoupl = Berendsen
    pcoupltype = isotropic ; uniform scaling of box vectors
    tau_p = 0.5 ; time constant (ps)
    ref_p = 1.0 ; reference pressure (bar)
    compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)

    ; VELOCITY GENERATION
    ;----------------------------------------------------
    gen_vel = yes ; Velocity generation is on (if gen_vel is 'yes', continuation should be 'no')
    gen_seed = -1 ; Use random seed
    gen_temp = 298.15

    ; FREE ENERGY
    ;----------------------------------------------------
    free-energy = yes
    sc-alpha = 0.5
    sc-power = 1
    sc-sigma = 0.3
    sc-coul = no
    init-lambda-state = 0
    coul-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
    vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
    nstdhdl = 100
    dhdl-print-energy = yes
    calc-lambda-neighbors = -1
  • 成品模拟mdp

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92

    ;====================================================
    ; Production MD
    ;====================================================

    ; RUN CONTROL
    ;----------------------------------------------------
    integrator = sd ; stochastic leap-frog integrator
    nsteps = 5000000 ; 2 * 5,000,000 fs = 10 ns
    dt = 0.002 ; 2 fs
    comm-mode = Linear ; remove center of mass translation
    nstcomm = 100 ; frequency for center of mass motion removal

    ; OUTPUT CONTROL
    ;----------------------------------------------------
    nstxout = 50000 ; save coordinates to .trr every 100 ps
    nstvout = 0 ; don't save velocities to .trr
    nstfout = 0 ; don't save forces to .trr

    nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps
    xtc-precision = 1000 ; precision with which to write to the compressed trajectory file
    nstlog = 1000 ; update log file every 2 ps
    nstenergy = 1000 ; save energies every 2 ps
    nstcalcenergy = 100 ; calculate energies every 100 steps

    ; BONDS
    ;----------------------------------------------------
    constraint_algorithm = lincs ; holonomic constraints
    constraints = h-bonds ; hydrogens only are constrained
    lincs_iter = 1 ; accuracy of LINCS (1 is default)
    lincs_order = 6 ; also related to accuracy (4 is default)
    lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
    continuation = no ; formerly known as 'unconstrained-start' - useful for exact continuations and reruns

    ; NEIGHBOR SEARCHING
    ;----------------------------------------------------
    cutoff-scheme = group ; 'group' is default in 4.6 - from GMX5 on 'Verlet' is default
    ns-type = grid ; search neighboring grid cells
    nstlist = 10 ; 20 fs (default is 10)
    rlist = 1.2 ; short-range neighborlist cutoff (in nm)
    pbc = xyz ; 3D PBC

    ; ELECTROSTATICS & EWALD
    ;----------------------------------------------------
    coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
    rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
    ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
    pme-order = 6 ; interpolation order for PME (default is 4)
    fourierspacing = 0.10 ; grid spacing for FFT
    ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb
    optimize-fft = no ; don't calculate the optimal FFT plan for the grid at startup

    ; VAN DER WAALS
    ;----------------------------------------------------
    vdwtype = Switch
    rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
    rvdw-switch = 0.9
    DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure

    ; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
    ;----------------------------------------------------
    tc_grps = System ; 温度耦合组,系统
    tau_t = 1.0
    ref_t = 298.15

    ; PRESSURE COUPLING
    ;----------------------------------------------------
    pcoupl = Parrinello-Rahman
    pcoupltype = isotropic ; uniform scaling of box vectors
    tau_p = 2 ; time constant (ps)
    ref_p = 1.0 ; reference pressure (bar)
    compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)

    ; VELOCITY GENERATION
    ;----------------------------------------------------
    gen_vel = yes ; Velocity generation is on (if gen_vel is 'yes', continuation should be 'no')
    gen_seed = -1 ; Use random seed
    gen_temp = 298.15

    ; FREE ENERGY
    ;----------------------------------------------------
    free-energy = yes
    sc-alpha = 0.5
    sc-power = 1
    sc-sigma = 0.3
    sc-coul = no
    init-lambda-state = 0
    coul-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
    vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
    nstdhdl = 100
    dhdl-print-energy = yes
    calc-lambda-neighbors = -1

总共包含31个状态,即init-lambda-state的值为0 ~ 30

xvg文件

采用Python脚本提取自由能:
https://github.com/choderalab/pymbar
https://github.com/MobleyLab/alchemical-analysis
比如:

1
python alchemical_analysis.py -f dhdl.*.xvg -s 500

其他参考资料

[1] ALDEGHI M, BLUCK J P, BIGGIN P C. Absolute Alchemical Free Energy Calculations for Ligand Binding: A Beginner’s Guide [M]//GORE M, JAGTAP U B. Computational Drug Discovery and Design. New York, NY; Springer New York. 2018: 199-232. DOI: 10.1007/978-1-4939-7756-7_11