计算化学公社

标题: g_mmpbsa计算带电底物与酶蛋白的结合自由能为正相关问题请教 [打印本页]

作者
Author: uenh1998 时间: 2024-2-13 23:14
标题: g_mmpbsa计算带电底物与酶蛋白的结合自由能为正相关问题请教
   各位老师晚上好，我想就g_mmpbsa计算带电配体与蛋白的结合自由能发现的问题来请教一下大家。
   复合物为酶-底物的体系。其中配体带1个正电荷，其每个原子所带电荷由RESP2方法得到。利用gromacs做蛋白-配体的分子动力学模拟，力场为AMBER14SB+GAFF。大概是酶-底物体系的缘故，配体相对蛋白的RMSD曲线无法保持非常稳定的值，这里截取RMSD曲线较为平衡的那一段来计算配体-蛋白的结合自由能，方法是g_mmpbsa。结果如图。

   但是结果显示结合自由能为正，意味着底物-酶处于不利结合相互排斥的情况，在酶-底物这个体系中出现这个情况合理吗？计算一定是错的吗？我是看到社长的 http://bbs.keinsci.com/forum.php?mod=viewthread&tid=27434&highlight=%CD%D1%C0%EB 这个帖子里最后一句话“比如之前我看到有人研究的是酶，问怎么小分子在里面没法一直呆住。实际上要求小分子呆住本来也没意义，本身底物小分子偶然进去之后被催化完就走掉了，何故期望小分子能稳定呆在里面？” 感觉我算的结合自由能为正应该没问题，但是审稿人针对这一点提问，因此前来请教。

   但之后从计算RESP2电荷步骤开始，到计算结合自由能的步骤重复了多次，都是这个情况，意味着我的计算应当没有不合理之处。那么针对审稿人提出的结合自由能为正的质疑，我该如何准备相关的材料从哪些角度去回答呢。

(, 下载次数 Times of downloads: 23)

作者
Author: sobereva 时间: 2024-2-14 11:39
尝试不同的溶剂模型（GB的不同有效Born半径计算方法、PBSA），也可以尝试gmx_mmpbsa (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/)算算看

作者
Author: uenh1998 时间: 2024-2-14 13:11
本帖最后由 uenh1998 于 2024-2-15 00:09 编辑

sobereva 发表于 2024-2-14 11:39
尝试不同的溶剂模型（GB的不同有效Born半径计算方法、PBSA），也可以尝试gmx_mmpbsa (https://valdes-tresa ...

感谢社长回复。我还有一些问题，想请社长解答一下我的疑惑：
1、所以，结合能为正，也不一定是不合理的吧。
2、以下是我的pbsa.mdp里的参数设置，想请教一下社长，您说的“GB的不同有效Born半径计算方法、PBSA”这个，是在哪里进行修改呢。

;Polar calculation: "yes" or "no"
polar = yes

;=============
;PSIZE options
;=============
;Factor by which to expand molecular dimensions to get coarsegrid dimensions.
cfac = 1.5

;The desired fine mesh spacing (in A)
gridspace = 0.5

:Amount (in A) to add to molecular dimensions to get fine grid dimensions.
fadd = 5

;Maximum memory (in MB) available per-processor for a calculation.
gmemceil = 4000

;=============================================
;APBS kwywords for polar solvation calculation
;=============================================
;Charge of positive ions
pcharge = 1

;Radius of positive charged ions
prad = 0.95

;Concentration of positive charged ions
pconc          = 0.150

;Charge of negative ions
ncharge = -1

;Radius of negative charged ions
nrad = 1.81

;Concentration of negative charged ions
nconc = 0.150

;Solute dielectric constant
pdie = 2

;Solvent dielectric constant
sdie = 80

;Reference or vacuum dielectric constant
vdie = 1

;Solvent probe radius
srad = 1.4

;Method used to map biomolecular charges on grid. chgm = spl0 or spl2 or spl4
chgm          = spl4

;Model used to construct dielectric and ionic boundary. srfm = smol or spl2 or spl4
srfm          = smol

;Value for cubic spline window. Only used in case of srfm = spl2 or spl4.
swin = 0.30

;Numebr of grid point per A^2. Not used when (srad = 0.0) or (srfm = spl2 or spl4)
sdens = 10

;Temperature in K
temp = 300

;Type of boundary condition to solve PB equation. bcfl = zero or sdh or mdh or focus or map
bcfl = mdh

;Non-linear (npbe) or linear (lpbe) PB equation to solve
PBsolver = npbe

;========================================================
;APBS kwywords for Apolar/Non-polar solvation calculation
;========================================================
;Non-polar solvation calculation: "yes" or "no"
apolar = yes

;Repulsive contribution to Non-polar
;===SASA model ====

;Gamma (Surface Tension) kJ/(mol A^2)
gamma          = 0.0226778

;Probe radius for SASA (A)
sasrad       = 1.4

;Offset (c) kJ/mol
sasaconst    = 3.84982

;===SAV model===
;Pressure kJ/(mol A^3)
press          = 0

;Probe radius for SAV (A)
savrad       = 0

;Offset (c) kJ/mol
savconst       = 0

;Attractive contribution to Non-polar
;===WCA model ====
;using WCA method: "yes" or "no"
WCA          = no

;Probe radius for WCA
wcarad       = 1.20

;bulk solvent density in A^3
bconc = 0.033428

;displacment in A for surface area derivative calculation
dpos = 0.05

;Quadrature grid points per A for molecular surface or solvent accessible surface
APsdens = 20

;Quadrature grid spacing in A for volume integral calculations
grid          = 0.45 0.45 0.45

;Parameter to construct solvent related surface or volume
APsrfm       = sacc

;Cubic spline window in A for spline based surface definitions
APswin       = 0.3

;Temperature in K
APtemp       = 300

作者
Author: sobereva 时间: 2024-2-15 13:50

uenh1998 发表于 2024-2-14 13:11
感谢社长回复。我还有一些问题，想请社长解答一下我的疑惑：
1、所以，结合能为正，也不一定是不合理 ...

1 不能说原理上不合理，但这种情况需要更严格的验证，避免由于计算模型不合适导致
2 我不用g_mmpbsa，细节不清楚

作者
Author: uenh1998 时间: 2024-2-16 00:04

sobereva 发表于 2024-2-15 13:50
1 不能说原理上不合理，但这种情况需要更严格的验证，避免由于计算模型不合适导致
2 我不用g_mmpbsa，细 ...

感谢社长的回复。目前还有一些小的疑问想请教一下社长老师。
1、这个如果需要排除计算模型不合适导致的问题，应该从哪些角度进行呢，还请社长指导一下

。
2、目前我只想到了计算一下配体的范德华表面的静电势着色图以及结合位点附近的静电势，如果表现为不利的相互作用，那么即可验证计算得到的正结合能是合理的。也不知道这个想法是不是对的

作者
Author: sobereva 时间: 2024-2-16 20:47

uenh1998 发表于 2024-2-16 00:04
感谢社长的回复。目前还有一些小的疑问想请教一下社长老师。
1、这个如果需要排除计算模型 ...

原子电荷、力场参数、溶剂模型（模型的种类以及具体参数。以前我写过《谈谈分子模拟中的隐式溶剂模型与GB模型》http://sobereva.com/42，里面有一些简要讨论）

光检验这个说明不了什么问题

作者
Author: uenh1998 时间: 2024-2-18 13:03

sobereva 发表于 2024-2-16 20:47
原子电荷、力场参数、溶剂模型（模型的种类以及具体参数。以前我写过《谈谈分子模拟中的隐式溶剂模型与GB ...

   感谢社长的回复。在这里我考虑了一下社长您给的这几个方面，打算给审稿人的回复如下：
   1.力场参数：g_mmpbsa计算需要准确的电荷参数来描述蛋白质和配体之间的相互作用。分子动力学模拟所使用的蛋白力场为AMBER14SB。如果电荷参数的选择不当或者参数化方法存在问题，可能会导致计算结果的偏差。
   2.原子电荷：配体分子的原子电荷参数也是影响结合自由能计算结果的重要因素之一。如果原子电荷参数的准确性受到质疑，例如由于参数化方法的限制或者参数化数据的不足，就有可能导致计算结果的偏差。对于配体，我们采用RESP2思想来计算配体的原子电荷。且“Schauperl M, Nerenberg P S, Jang H, et al. Non-bonded force field model with advanced restrained electrostatic potential charges (RESP2)[J]. Communications chemistry, 2020, 3(1): 44.”文献说明了RESP2电荷思想比以往的RESP方法多数情况下还会更好。因此对于该配体，采用RESP2思想是没有问题的。
      3.溶剂模型：结合自由能的计算涉及溶剂模型的选择，包括溶剂的介电常数等参数。不同的溶剂模型可能会对计算结果产生不同的影响。g_mmpbsa计算通常使用隐式溶剂模型来描述溶液环境中的溶质溶解行为。虽然隐式溶剂模型具有计算效率高的优点，但它们也是对真实溶液环境的简化描述，可能无法准确地捕捉溶剂分子的动态行为和局部结构；隐式溶剂模型可能忽略了溶剂分子的极化行为，或者采用了静态的极化模型。在某些情况下，溶剂分子的动态极化行为可能对溶质的溶解行为产生重要影响，特别是在配体的结合界面附近。且隐式溶剂模型通常需要指定介电常数来描述溶液环境的极性程度。介电常数的选择可能会影响计算结果的准确性。而且，PB模型与显式水模型仍有一定差距，将水分子进行了“连续化”的近似，故不能表现与溶剂的强相互作用（如氢键）、特殊相互作用（如水桥）。由于PB模型下的近似并不抗衡离子，所以对于带电体系（如我们研究的带电配体），不能表达抗衡粒子与它的相互作用，这会高估静电相互作用能，从而导致最后的结合自由能偏大，甚至为正的结果。

社长觉得这样的回复有没有说服力呢；
此外我还在想需不需要其它的一些计算来支持这些回复，但是还没想到该计算一些什么作为补充。想问下社长有没有什么建议呢。感谢社长过年期间在回答我的帖子，社长新年快乐！

欢迎光临计算化学公社 (http://bbs.keinsci.com/)