隐式溶剂模型.docx
- 文档编号:27967826
- 上传时间:2023-07-06
- 格式:DOCX
- 页数:9
- 大小:23.03KB
隐式溶剂模型.docx
《隐式溶剂模型.docx》由会员分享,可在线阅读,更多相关《隐式溶剂模型.docx(9页珍藏版)》请在冰豆网上搜索。
隐式溶剂模型
1.Poisson-Boltzmann,PB模型============================
在量子化学计算和使用经典分子力学进行模拟中,常会用到隐式溶剂模型,特别是作为一种快速估算溶剂化自由能的方法,在蛋白-底物分子作用计算,和蛋白质折叠等方面特别有用。
量子化学gauss中,还会经常用来估算溶剂中某些分子的pKa。
那么,隐式溶剂如何做到不使用溶剂而估算溶剂的影响呢?
它的原理是什么?
我们首先从最基本的问题谈起,即溶剂化效应是怎么回事。
溶剂化效应是发生在溶剂中的很多化学现象的关键。
它影响反应的自由能势能面,改变反应速率和热力学平衡。
我们知道,一旦知道一个反应的势能面,我们就知道了它的全部热力学和动力学信息。
举个例子,比如乙烷的二面角的旋转,如果我们把二面角作为反应坐标,则如果我们知道了这个旋转势能面,那么我们就知道了每个二面角取值的能量(势能),那么根据Boltzmann分布,我们可以计算任意温度下,乙烷对各二面角取值的分布概率:
{P_{torsion}}=exp(-E_{torsion}/kT)/P_0}{P_{torsion}}是能量为{E_{torsion}}的构象的权重,而P0是配分函数,P0是把所有构象的权重加起来。
这个概率就是所谓的热力学平衡分布,即体系平衡后某构象的浓度等于总浓度乘以该构象的概率:
{C_{torsion}=C_{total}*P_{torsion}}而动力学,就是翻越构象能垒的速率,可由能垒高度ΔE得到:
反应速率常数:
{k=exp^{-\DeltaE/kT}}反应速率由阿伦尼乌斯公式,即产物浓度变化:
{d/dt=A*k*C_{reactant}}(C_reactant为反应物浓度)由上面可知反应势能面真的很厉害。
有了反应势能面,别无他求。
有些初学者往往觉得跑个Gauss没什么,做MD才给力。
其实,如果量子化学计算能够计算势能面,干嘛要跑MD呢?
MD不过是模拟体系在那个势能面上随时间的变化而已。
不过,量子化学计算个单分子在气相中的反应,搜个过渡态啥的,很准确。
但是对于溶剂中的反应呢?
这时候,溶剂的作用,即溶剂分子在反应物分子周围来来往往的影响就不能不考虑了。
回到上面乙烷的例子,如果在水里面旋转,肯定要受到水分子的阻力什么的影响。
那么多的水分子,随便哪个分子动一下,阻力就不一样,这就是熵效应,即所有水分子的不同构象对乙烷旋转势能面的影响。
在考虑了周围溶剂的影响后,乙烷的势能面已经不叫势能面了,而叫旋转自由能势能面。
自由能包括内能和熵两个部分。
当然,熵不仅包括溶剂的,也包括乙烷自身的构象熵,比如它自身在水里面的振动,转动,平动,都影响最后的自由能势能面结果。
那么,怎么考虑溶剂化效应对反应的势能面的影响呢?
即,如果我们要计算自由能势能面,那有2个方法:
第一,跑MD,随机模拟溶液体系的一些构象,然后根据体系中乙烷取各构象的概率,反推出自由能势能面;第二,不跑MD,还是用量子化学的方法来计算单个乙烷的旋转势能面,但是使用隐式溶剂模型,即使用一个经验的公式,来把溶剂的影响的平均结果考虑进来,这样,直接得出自由能势能面。
怎么平均化溶剂作用影响呢?
学过高中物理的都知道,设置一个等效介电常数。
即,真空中两电荷作用力{F=(1/4\pi\epsilon_0)*q_1*q_2/r^2}
那,放到溶剂中,{F=1/4\pi\epsilon0*\epsilon_r*q_1*q_2/r^2}比如,水的{\epsilon_r=78}当然,我们在拉格朗日图景中,要使用能量而不是牛顿力学中的力。
即,我们要求的双点电荷的Coulomb势能U(r12):
{U(r_{12})=(1/4\pi\epsilon_0)*q_1*q_2/r}十分简单。
But,waitaminute...实际的隐式溶剂模型,并不使用Coulomb公式,而是使用Poisson-Boltzmann公式!
!
!
为什么呢?
模型不一样。
Coulomb公式用的是点电荷模型,q1,q2,距离r。
但是,对于真实的分子,是由原子核加核外电子云构成,电子云是有体积的。
量子化学计算就是核外电子云的密度变化。
即使是使用经典分子力学力场,AMBER,CHARMM什么的,也是把原子看成是带电荷的,有半径r的球的,这样,必然每个原子附近都有不同的电荷密度。
如果把原子看作点电荷,进行对加和计算,则收敛慢且截断误差很大。
一个更聪明,更快速的办法是:
把体系分成均匀格点,计算每个格点的静电势Φ(r),然后在该点上的电荷q(r)所具有的能量则为{U(r)=q(r)*\phi(r)}这样,最后把各点上的U(r)加起来,就得到总的静电势能U(注意静电势与静电势能的差别)。
这就是Poisson方法的思想。
Poisson与Coulomb对加和计算的差别,类似于量子化学计算中DFT方法与HF的差别,PB把Coulomb的多体作用转化为随空间变化的密度问题,势能等于密度乘以静电势。
现在,最关键的,是静电势如何得到?
它是由Poisson解决的。
简单说,Poisson给出了一个微分方程,这个微分方程描述了介电常数,电荷密度分布和静电势的关系,给出了一个由一大片电子云或几个空间分布的带电原子球ρ(r)计算静电势Φ(r)的关系。
Poisson公式:
{\nabla=\rho(r)}?
:
拉普拉斯算符,可理解为一阶偏导r:
表示空间格点,你可以把空间用线切割开,把r想象成一个个格点位置(注意跟上面的表示两电荷距离的r定义不同),格点是使用有限差分方法求解微分方程的常用方法ε:
介电常数Φ:
所要求的静电势ρ:
各格点上的电荷密度分布这里,Φ(r)是代表体系r处的静电势,是我们要求的量,它的二阶导由该点的电荷密度ρ(r)及其介电常数ε(r)决定。
虽然这个公式的物理意义我理解不了,但是我只需要知道通过它我们可以得到静电势Φ(r)就足够了,即由上式可得:
{\nabla^2\phi(r)=-4\pi*\rho(r)/\epsilon(r)}
这样,我们得到在各点r处静电势能的二阶导?
^2Φ(r),最后对整个盒子所有格点进行二重数值积分即可得到最终静电势,乘以各点上的电荷得到静电势能U(r)=q(r)*Φ(r)。
MD中一般使用的计算静电势能的PME方法,即ParticleMeshEwald方法,原理就是利用Poisson公式进行格点积分。
上面只说了Poisson公式,不是说隐式溶剂用的是Poisson-Boltamann公式吗?
Boltzmann干了啥?
问题是这样的,上面的模型,用来计算蛋白质在一定浓度电解质溶液中的静电势能,存在一个系统误差,原因是没有考虑反离子的分布的贡献。
蛋白质往往是带电荷的,结果在真实溶液体系中,在一个蛋白的不同表面处,往往分布着较高浓
度的相反电荷。
这叫离子氛模型。
它造成的结果是,溶液中的静电作用的衰减要比由库仑公式计算的势能面的衰减快得多。
这叫离子屏蔽效应。
电解质溶液理论中使用Debye长度来屏蔽效应的强弱。
但是Debye长度是一个经验量,无法很方便地定义。
怎么能考虑这个屏蔽效应呢?
一个直接的想法是,把溶质比如蛋白附近的反离子浓度计算出来。
即一个改变了的电荷分布:
ρ(r)+Δρ(r),其中后面一项Δρ(r)就是溶液体系中反离子的浓度在r处的分布。
即:
{\nabla=\rho(r)+\Delta\rho(r)}通过这个改进的Poisson公式,可以更好的考虑蛋白质在电解质溶液即离子溶液中的静电势能。
当然,由于反离子分布本身改变静电势能,而静电势能又改变反离子浓度分布,因此,这必然是个自洽计算。
那么怎么计算反离子浓度呢?
Poisson-Boltzmann公式中的Boltzmann部分:
{\Delta\rho(r)=n_i(r)=n_i^0*exp(-q_i*\phi/kT)}{n_i(r)}:
溶液中离子种类i在r处的数密度(可以是分数个){n_i^0(r)}:
溶液中离子种类i的平均数密度{q_i}:
该种类离子的电荷{\phi(r)}:
盒子空间中该区域的静电势上式的物理意义是,某离子在盒子内某处数密度对平均密度的偏差,与其在该处的造成的静电势能的大小成自然幂反比关系。
换句话说,在蛋白附近处静电势Φ越高,则该处的反离子浓度越高,同号离子的浓度越低。
因为前者造成能量下降,后者造成能量上升。
对于一般尺寸的盒子,由于在数千到数万个格点上进行含指数计算的自洽计算,收敛相当相当慢,伤不起,只好对该指数进行tylor展开为多项式,如果截取多项式到第一项,则ni与Φ为一个线性关系,这样的近似称为线性PB计算LPBE.LPBE对于不带净电荷的小分子来说是足够准确的;但是,对于带电荷的,特别是电荷浓度较大的体系,通常要取到第二项,ni是Φ的三次方的函数,这样的近似称为非线性PB计算NPBE。
这方面的论文很多,尝试用各种妖娆的数值方式加快计算。
简单了解一下上面这些名词很有必要。
在很多采用一阶近似的线性PB模型中,为了避免屏蔽效应考虑的不充分,会在Boltzmann那一部分的前面加上一个系数{K'^2=8\pi*N_A*e^2*I/1000kT},直接考虑Debye-Hukel效应。
其中{Debyelength=1/k=1/\sqrt{(K'^2/\epsilon)}}
2.GernalizedBorn-GB模型===========================
为了计算一个半径为r,带净电荷为q的球的溶剂化自由能,即把它从气相移到介电常数为ε的溶剂中,Born提出了一个及其简单的计算公式:
{\DeltaG_{elec}=-(q^2/2r)*(1/\epsilon-1)}可以简单理解为,如果你想往溶剂里移动了一个电荷q,那么溶剂必然会被极化,在电荷所在位置产生了一个同号的q,你的移动必须克服这电荷与镜像之间的排斥力。
这个模型的奇妙之处在于只有一个可调参数,即半径r.一开始born使用原子的共价半径r,计算结果虽然误差较大,但定性准确。
考虑到模型的简单性,可以说是出乎意料。
后来有人用vdW半径计算,误差大大减小,但仍需提高。
于是,不断有人试图使用不同的半径,经验关联的,量子化学计算的。
总而言之,Born
半径是Born模型的关键。
GeneralizeBorn模型,即广义Born模型,或一般化了的Born模型,是用它来估算分子的溶剂化自由能,即把分子中的每个组成原子用上面的模型估算一下,然后加和得到总的溶剂化自由能;你会问,分子表面的原子与溶剂的作用比较强,分子内部的原子与溶剂几乎没有作用,怎么能加和呢?
不要紧,给分子不同位置的原子设置不同的有效半径即可,这个有效半径可以用来反应原子离表面的距离,或者说被埋的程度。
另外,你必须注意到,每个原子不仅与自身的镜像作用,而且还会与其它原子作用。
{\Delta
G_{elec}=-\sum\{(q_i^2/f_i)*(1/\epsilon-1)\}-\sum\{(q_i*q_j/f_{ij})*(1/\epsilon-1)}\}
前面一项就是原来的Born,后面一项是原子间两两作用。
这两项中的r都被f取代了。
f的计算:
{f=\sqrt{r_{ij}^2+a_{ij}^2*exp^{-D}}}f等于原子间距rij和原子born半径aij按指数衰减的均方根;{a_{ij}=\sqrt{a_i*a_j}}其中有效半径aij等于两原子的有效半径的几何平均;{D=r_{ij}^2/(2a_{ij})^2}D是f中aij^2按指数衰减的幂,这个衰减类似于考虑周围反离子分布的屏蔽效应;GB模型是PB模型的特例。
(Tobeaddressed)当然,这里仅仅计算了静电作用对溶剂化自由能的贡献。
还有一部分贡献来自于分子体积破坏了水分子间的氢键网络,从而产生的能量代价。
这一部分的贡献对于中性分子的溶剂化自由能大概占10%。
怎么计算这一部分的能量呢?
先问破坏了H键网络的后果是什么呢?
水分子在分子表面形成了表面张力,由朗格谬而,表面形成能等于表面张力乘以表面积,所以这部分自由能必然是表面积的线性函数。
所以,要计算的就是溶剂可及表面solventaccessiblearea,SA。
由于一般使用vdw半径作为探测溶剂可及表面的基准,所以可称之为ΔGvdw。
即:
ΔGvdw=σ*A最终:
{\DeltaG_{solvation}=\DeltaG_{elec}+\DeltaG_{vdw}}因此,这种估算溶剂化自由能的方法称为GB/SA模型,分子模拟中有时称为MM/GBSA。
必须注意有报道称MM/GBSA模拟多肽的构象得到的结果与显式溶剂模型的结果有显著的不同(参考英文wiki:
implicitsolvation-GBSA条目),其MM/GBSA不能正确识别蛋白native的结构。
还有盐桥过于稳定的毛病,可能是因为屏蔽效应估计不足;alpha-helix的分布也明显高估。
看来,这个PB的近似模型尽管简单,但是要慎重使用阿。
要想提高计算精度,可以使用PB/SA模型,即静电部分的贡献,占溶剂化自由能的主要部分(>90%),使用PB来计算。
当然,这一部分要在几千个格点上做自洽计算,所以要慢的多。
参考文献
1.Wiki:
implicitsolvation
2.D.R.Livesay,ContinuumElectrostaticMethods,~drlivesa/ITSC8311/ITSC8311_elec.pdf
3.P.Grochowski,J.Tryska,Continuummolecularelectrostatics,salteffects,andcounterionbinding--areviewofthePoisson-Boltzmanntheory
anditsmodifications.Biopolymers.(2008)89,93-113.
本文主要讨论隐式溶剂模型和GB模型(广义Born模型)的各个方面。
隐式水模型优势主要有下:
减少体系的粒子数从而降低计算开销。
没有溶剂摩擦效应,可以加快构象变化。
不显著、具体地表现溶剂运动而直接表现平均效果,便于计算自由能,不需要动力学很长时间采样。
不需要显式水的预平衡过程。
适用于常PH模拟、REMD等方法。
适合计算能量面,结果比较光滑,真实的水会带来噪音,溶剂的运动、碰撞使主体分子产生许多能量局部极小点。
等等......
总体来说,隐式水模型相对于显式水模型,是一个离散化到连续化的过程。
溶剂模型的近似关系如下(越左边越真实,越右边近似程度越大):
显式溶剂模型-->隐式溶剂模型-->明确拆分为静电、非极性贡献-->PB-->GB
分子在溶剂环境下的总能量Etot=溶剂化能deltaG(solv)+溶质在真空条件下的能量Evac。
(对于极化力场,溶剂对溶质的极化作用造成的能变也算进在deltaG(solv)里。
)
deltaG(solv)=deltaG(elec)+deltaG(nonpolar)=deltaG(elec)+deltaG(vdw)+deltaG(cavity)
注:
有时明确包含氢键项。
其中deltaG(elec)就是溶剂环境下的静电作用能与真空下的静电作用能之差,主要包括溶剂与溶质间的静电作用,以及溶剂对溶质原子间静电作用的屏蔽效果。
其中deltaG(nonpolar)就是除此以外,即体系各原子皆无电荷时,溶剂环境与真空环境下能量之差。
它又分为deltaG(vdw)和deltaG(cavity)。
deltaG(vdw)是溶质与溶剂的范德华作用,由于范德华作用随距离衰减快,所以只考虑第一溶剂化层的溶剂,此项造成能量降低。
deltaG(cavity)表现的是溶质反抗溶剂压力形成洞穴的做功以及溶质周围溶剂(第一溶剂化层为主)的重构造成的熵罚,此项造成能量升高。
因为deltaG(nonpolar)的两项都主要涉及第一溶剂化层,它正比于SASA,故最简单、常用的表达是a*SASA+b,而比例系数a由拟合实验值得到,b一般为0。
另有其它方案计算,比如包括溶质体积、溶质直径等信息在内的方程。
也有比如
deltaG(vdw)通过SAS的积分得到,而deltaG(cavity)仍用正比于SASA的方案。
在GB模型下的MD模拟中,对于原始状态的生物大分子,由于SASA变化很小,原子因非极性作用导致的受力往往可以忽略,若需计算的话则是求它对坐标的导数。
由于a*SASA+b的a为正值(一般为7.2cal/mol/?
^2),可见deltaG(nonpolar)的效果是减小溶质的SASA,此项体现了疏水效应。
deltaG(elec)项一般通过GB或PB方法计算,GB速度快但一般准确度低于PB。
目前常用的GB模型计算方程由Still于1990年提出(原文),依据是库仑定律和born方程,推导出另一种形式。
当粒子相离较远时,方程等价于库仑定律+born方程能量,而相距为0时等价于born方程,而当距离较近时等价于Onsager反应场能量。
方程中还可以再引入单价粒子的屏蔽校正项。
GB方程中重要参数是有效born半径。
某个原子电荷与溶剂作用对deltaG(elec)的贡献实际上是把分子中其它原子电荷去掉后,溶剂化状态与真空状态能量差。
如果用的有效born半径以GB方程算的结果与之一样,就是正确的born半径。
或者粗略来说:
某原子用有效born半径算得的GB式中的“自身项”的能量应当与真实的分子中此原子与溶剂的静电相互作用能一致。
这样,有了正确的born半径,就可以用GB式得到合理的deltaG(elec)。
有效born半径对计算deltaG(elec)有重要影响,模拟过程中分子构象不断改变,依赖于它的有效born半径也不断改变,需要每一步重新计算,若构象变化不明显可每隔一定步数重新计算,这是GB计算量中的主要部分之一。
计算有效born半径在不同程序中有不同方法。
一般使用CFA(库仑场近似),有效born半径就可以对溶质分子表面内的空间进行积分获得。
而分子内空间的范围不便于描述,简单的方法是直接将原子VDW球内区域叠加作为分子内空间,称GB-HCT,但这样就忽略了原子VDW空间之间的缝隙,可以乘以4/3来近似校正。
也有人根据原子被埋程度的概念提出了基于经验的简单、快速的函数,其中引入可调参数,在amber里称为GB-OBC方法,结果明显优于GB-HCT,应注意可调参数是通过预先优化得到的,面向不同体系模拟应使用不同参数。
amber中还支持更新的GBn方法,结果优于GB-OBC,适合蛋白而不适于核酸体系。
GB方程的优点是形式简单,计算快速,结果较准确,而且有简单的解析导数,可直接计算GB环境下原子在MD过程中的受力。
GB可以结合分子力场及量子力学模拟各种分子体系的溶剂化效果。
可以用在类似MM/PBSA的结合自由能计算的静电项中,即MM/GBSA。
一些模拟方法则只能使用GB,比如amber中可以实现的结合MC的常PH模拟,原理是根据Metropolis判据在MD过程中动态决定是否某点被质子化或去质子化。
对于REMD,因为随着体系粒子数的增加,需要的副本数目飙升(否则能量重叠较差,交换概率低而起不到效果),所以总是结合GB方法来显著降低粒子数量。
对于不大的体系,GB比显式溶剂模型有更快的速度。
但GB方法也有缺点。
与显式水模型仍有一定差距,将水分子进行了“连续化”
的近似,故不能表现与溶剂的强相互作用、特殊相互作用(如水桥)。
对于膜体系,溶剂、溶质、膜都有着不同的介电常数,这样复杂的情况介电常数环境下GB也难以表达。
对于电荷较多的核酸体系,也不能较好表达多价抗衡粒子与它的相互作用。
一些研究也表明用GB模拟多肽会产生与显式溶剂模型不同的构象,与实验结果有一定偏差。
计算deltaG(elec)项时,相对于原理严格的PB也有差距(即便使用最佳的有效born半径),但是好处是GB计算速度快得多。
不过PB的结果亦不可作为绝对的金标准,因为一些误差在PB引入的近似中就已经出现了。
在计算速度上,显式溶剂体系计算量是O(N^2),但使用Ewald等方法后,可降至Nlog(N)以下,而GB模型如果在cutoff上不做近似(往往取无限大),对于大体系反倒更慢。
除GB、PB以外也有其它一些方案,最简单的隐式模型是介电常数依赖于作用距离的方法,其中库仑作用的介电常数等于粒子间距离(即静电作用能1/r变为了1/(r^2)),以体现溶剂的屏蔽效果,此方法精度显然低于GB。
近来还有一些新方法,如AGBNP(AnalyticalGeneralizedBornplusNonpolar),在GB基础上引入另外形式的deltaG(nonpolar)项。
ALPB(analyticallinearizedPoisson-Boltzmann),其方程类似GB(在无限介电常数下回归为GB方程),故有着基本一致的计算速度,但引入了有效分子静电尺寸参数,模拟水比GB有着更好的精度,结果更接近于PB。
这些新模型理论上有着更好的结果,但仍需广泛地应用于模拟来检验。
还有隐式与显式溶剂模型相结合的方法,也就是在主体分子外面包一层溶剂分子,往往以弹簧势限制住溶剂避免跑走,而更外面则用隐式溶剂模型,结合了两种方法的一些优点。
目前Amber对隐式溶剂模型支持较好,支持GB、ALPB,也有内建的PB模块pbsa(老版本挂的是第三方的delphi)。
而Gromacs在最新版本4.0.5中尚不支持GB/PB,需要外接第三方软件如APBS,但在接下来的版本中即将加入GB。
对于GB模型,适用领域与不适用领域的界限尚为模糊,一般来说,对于必须用GB的地方,或者已证明有明显优势的地方可以使用,但如果不确定,尽量还是用显式溶剂模型,毕竟GB是基于多层近似后的溶剂模型。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 溶剂 模型