1 引言
水分子作为最常见的溶剂与其他众多化合物相互作用,在分子水溶液的研究和分子力场开发中处于核心位置。水分子体系不论是在量子力学模拟中,还是在经典分子模拟中,是研究的热点之一。水是细胞内的主要溶剂(占细胞体积70%~80%),生物分子(如蛋白质、核酸)在水溶液中才能维持天然构象与功能活性(如酶催化、DNA复制)。非水环境(如真空)会导致蛋白质异常折叠或功能丧失。因此,水是模拟生物分子在生理环境中动态行为的重要溶剂,选择合适的水模型可以提高模拟的精度。
分子动力学模拟是研究蛋白质、核酸等生物大分子性质不可缺少的工具
[1 -3 ] 。然而,它们的有效性,取决于描述所有原子之间作用力的分子力场的准确性。因此,发展新的力场用于研究生物分子水溶液体系的结构和变化,是十分必要的。基于分子力学的基本近似,只要某一势能函数能够很好地描述体系内的各种相互作用,体系内的原子和分子的运动,就可以用源于该势能函数的经典力学来描述,该势能函数以及相关的参数就称为分子力场。早期传统分子力场(FFs)为了能以适度的计算能力模拟生物分子体系,通常使用原子上的固定电荷作为参数,这导致固定电荷力场是有局限性的
[4 -9 ] ,因为它忽略了带电原子之间的极化和电荷转移,而这两点对生物系统的识别和稳定性至关重要
[5 ,9 ] 。与固定电荷力场相比,极化力场中的电荷动态变化,能够更好地体现极化效应。在过去的几十年中已经开发出多种可极化力场
[10 ] ,包括AMBER
[11 -13 ] 、AMOEBA
[14 -15 ] 、CHARMM Drude
[6 ] 、CHARMM浮动电荷
[16 -20 ] 、SIBFA
[21 ] 、GEM
[21 ] 和ABEEM
[22 -31 ] 。因此,发展极化力场用于生物分子体系的分子动力学模拟是十分必要的。
本综述总结了ABEEM极化分子力场(ABEEM PFF)在生物分子水溶液的分子动力学模拟方面的相关工作。介绍了ABEEM-7P水模型及ABEEM PFF在生物分子体系中的应用,通过一些分子动力学模拟实例,证明ABEEM PFF可以准确预测生物分子在水溶液中的结构、电荷分布、能量等动态性质。最后对ABEEM PFF未来的发展进行了展望。
2 理论方法
2.1 定态Schrödinger方程的Born-Oppenheimer近似求解
分子常是含有若干原子核和多个电子的体系,其波函数满足定态Schrödinger方程,该方程是分子中电子结构的基本方程。依据当分子构型有一微小变化时迅速运动的电子能立即进行调整的思想,Born和Oppenheimer简化分子体系的定态Schrödinger方程,把分子中核运动和电子运动分离,称为定核近似或B-O近似。在这个近似下,分子体系的波函数可以写为电子运动波函数 Ψ e R , r 与原子核运动波函数 Ψ n R 的乘积形式,与之相对应,分子体系的非相对论的定态Schrödinger方程可分解为两个方程,是量子化学电子结构理论中最重要的方程。分别是关于电子运动的式(1) 和关于原子核运动的式(2) 。
∑ i ( - 1 2 ∇ i 2 ) + V e R , r Ψ e R , r = E e ( R ) Ψ e R , r (1)
∑ A ( - 1 2 M A ∇ A 2 ) + ∑ A ∑ B < A Z A Z B R A - R B + E e R Ψ n R = E Ψ n R (2)
式中: R A 代表第A 个原子核的坐标; r 代表所有电子的坐标; R 代表所有原子核的坐标; Ψ e R , r 代表在某个核坐标 R 下的电子波函数; Ψ n R 代表原子核的波函数,E 代表体系的总能量。显然,在一定原子核坐标 R 下,体系电子的总能量 E e R ,由式(1) 决定,同时也出现在式(2) 中,称作势能面。在讨论分子结构、化学反应等问题时,势能面都是不可或缺的。
2.2 分子力学(MM)
MM是一种建立在经典力学基础上的计算方法。该方法使用经验或半经验参数计算分子的结构和能量,通常被称为力场(Force field,FF)方法。与量子力学方法不同,它不需要求解体系的波函数,而是基于Born-Oppenheimer近似,将电子和原子核的运动分离开来,只考虑原子的运动,并将体系的能量视为原子核位置的函数。分子力场描述的势能面,一般包括键长的伸缩振动势能、键角的弯曲振动势能、二面角的扭转势能、非共面二面角势能以及静电和范德华相互作用等部分。这些力场参数通常通过实验数据或量子化学计算拟合获得。使用力场方法可以计算大分子体系的多种性质,包括能量、结构和热力学性质等。这种方法为研究分子性质提供了有力的工具,特别是在探索大分子体系的行为时,具有广泛的应用。
2.3 常见力场
20世纪70年代,计算机模拟开始广泛应用于有机和生物分子体系的研究。分子动力学和蒙特卡罗统计力学开始用于模拟液态水、含有简单分子和离子的水溶液、有机液体以及真空中的小蛋白质等体系。这些工作中重点在于对分子内和分子间能量的准确描述。鉴于体系的大小和可用的计算资源,通常采用经典力场进行模拟,这些力场是有机分子的分子力学的研究基础,不同力场之间的差异主要在于所采用的力场模型的不同,包括是否考虑氢键相互作用或极化效应等因素。此外,各种力场还存在不同的数学函数公式和参数设置,这些因素导致了力场的多样性。力场在模拟分子相互作用时根据处理静电相互作用的不同可分为两种主要类型:固定电荷力场和极化力场。固定电荷力场的特点是在动力学模拟过程中,分子体系中各原子的电荷是固定不变的。这意味着即使分子的结构发生变化,原子上的电荷分布依然保持不变。由于它未能准确考虑分子中电荷的极化和分布变化,这种模型在模拟分子时存在一些局限。当涉及高电量离子之间的相互作用、离子与π
- 电子体系的相互作用以及离子的导电性和渗透性时
[14 ] ,使用极化力场能够更好地描述这些情况
[32 ] 。
2.3.1 固定电荷力场
固定电荷力场(Fixed-charge force field,FFF)概念的提出可以追溯到20世纪60年代末和70年代初,其核心思想是将分子中的原子电荷视为定值,这些电荷值的获取主要来源于QM计算和实验数据。AMBER(Assisted Model Building with Energy Refinement)力场
[33 -34 ] 作为最早的固定电荷力场模型之一,主要用于模拟蛋白质和核酸体系。随后,CHARMM(Chemistry at HARvard Macromolecular Mechanics)
[35 -36 ] 、OPLS(Optimized Potentials for Liquid Simulations)
[37 -39 ] 和GROMOS(GROningen MOlecular Simulation)
[40 -41 ] 等力场模型相继出现。
AMBER力场是由Kollman等
[33 ] 开发的力场,其力场的优势在于能快速并较准确地计算生物大分子如包含蛋白质、核酸和多糖等体系的性质,且力场参数较少。在AMBER力场的势函数体系中,总势能可分解为成键相互作用和非键相互作用,其展开式如
式(3) 。
E = ∑ K 2 ( b - b 0 ) 2 + ∑ K θ ( θ - θ 0 ) 2 + ∑ 1 2 V 0 1 + c o s ( n ϕ - ϕ 0 ) + ∑ ε ( r * r ) 12 - 2 ( r * r ) 6 + ∑ q i q j ε i j r i j + ∑ C i j r i j 12 - D i j r i j 10 (3)
由Karplus课题组
[35 ] 开发的CHARMM力场在处理于有机小分子和溶液中的大分子体系时能够获得与实验结果几乎一致的分子结构、能量等性质。其力场参数来自于实验值和从头算方法,CHARMM力场的势函数如
式(4) 所示,每一项分别代表键长势、键角势、Urey-Bradley相互作用、二面角势、非共面作用势、色散作用、范德华作用和库仑相互作用。
E = ∑ k b ( r - r e q ) 2 + ∑ k θ ( θ - θ e q ) 2 + ∑ k U B ( S - S ) 2 + e q ∑ V ω 1 + c o s ( n ω - ω e q ) + ∑ k χ χ - χ e q 2 + ∑ ε i j ( R i j r i j ) 12 - 2 ( R i j r i j ) 6 + ∑ q i q j 4 π ε i j r i j (4)
为了更好地模拟液体的热力学性质,Jorgensen等
[37 ] 开发了OPLS力场。该力场包含全原子和联合原子力场两个版本,主要针对模拟含溶剂体系,其成键项参数取自AMBER力场,非键项的参数由实验进行标定。
目前,除AMBER、CHARMM和OPLS等力场外,GROMOS力场也是用于生物分子模拟的最重要的固定电荷分子力场之一。该力场是通过精心设计而提高了计算效率的联合原子分子力场,适用于烷烃、蛋白、核酸凝聚相的模拟,势能函数形式如式(5) 。
E = ∑ 1 2 k l ( l 2 - l 0 2 ) 2 + ∑ k θ ( c o s θ - c o s θ 0 ) 2 + ∑ 1 2 k m 1 + c o s δ c o s ( m ω ) 2 + ∑ 1 2 k i ( ξ - ξ 0 ) 2 + ∑ 4 ε i j ( σ i j r i j ) 12 - ( σ i j r i j ) 6 + ∑ q i q j 4 π ε 0 r i j (5)
计算简单是固定电荷力场最突出的优点,使其适用于大规模的分子模拟,并在分子动力学模拟中应用广泛。例如,在研究药物设计时,固定电荷力场最突出的贡献在于可以预测药物分子与蛋白之间结合的靶点,从而可以对潜在的药物分子进行筛选。在研究蛋白质折叠和酶催化时,固定电荷力场可以从分子水平层面提供重要的计算结果。此外,固定电荷力场对于材料化学中界面现象和纳米材料的相变行为和热力学等性质也有研究。
在不断的改进和发展中,这些力场模型逐渐展现出各自独特的风格特点和应用。尽管固定电荷力场在许多领域中表现出较好的应用,但仍存在一些不足之处和局限性。由于固定电荷力场中原子电荷固定不变,并不考虑电荷的极化效应。因此在描述某些情况下分子间的相互作用时可能并不准确。例如,对极化效应的忽略使得气相下计算的中心分子电荷分布与跟溶剂计算得到的结果没有差异,导致构象的预测结果错误。当处于强电场或高极性环境下,模拟结果可能会偏离实际情况。为了克服这些局限性和更准确地描述分子间的相互作用,更为复杂的可极化力场被开发。
2.3.2 极化力场
极化力场(Polarizable force fields, PFFs)能够在模拟过程中显式处理实际环境中产生的极化效应。极化力场主要通过两种模型表示极化效应:第一种模型通过引入诱导偶极(多极子)
[15 ] 或Drude振子(也称为弹簧上的电荷模型或者壳模型)
[6 ,42 -43 ] 体现每个原子内的电荷重新分布;另一种是基于电荷在原子间转移的浮动电荷模型
[17 ,25 ] 。在诱导偶极模型中,每个极化位点在外电场
E 的作用下产生诱导偶极矩
μ ind ,
μ ind =
αE ,其中
α 是极化率
[15 ] 。这些诱导偶极也是总电场的一部分,因此极化不是加和性的,通常通过自洽迭代的方法进行确定。在Drude振子模型中
[6 ,42 -43 ] ,通过简谐弹簧将携带部分电荷的Drude粒子连接到原子上。Drude粒子的位移会形成诱导偶极矩。尽管在MD模拟过程中已采用扩展的Lagrangian方法来逼近精确解,但Drude粒子的位置也需要通过自洽迭代进行求解,以确保Drude粒子处于基态。浮动电荷模型是基于电负性均衡原理,通过每个位点电负性或化学势均衡来求得原子电荷。诱导偶极和Drude振子模型中允许在每个位点产生一个偶极。由于通常情况下Drude振子的位移非常小,故将其近似为一个点偶极。浮动电荷模型允许原子间电荷转移。简单的浮动电荷模型未考虑平面外的极化或单原子离子的极化。这些效应在浮动电荷模型中可通过在每个原子周围引入额外的虚位点表示
[28 ] 。
目前,应用于生物体系的极化力场包括AMBER
[11 -12 ] 、AMOEBA
[14 -15 ] 、CHARMM Drude
[6 ] 、CHARMM-FQ
[16 -17 ] 、SIBFA(Sum of interaction between fragments
ab initio )
[21 ] 、GEM(Gaussian electrostatic model)
[21 ] 以及ABEEM PFF
[22 -23 ] 等。
AMBER ff02pol
[12 ] 是最早用于蛋白和核酸的极化力场之一,其采用点电荷和简单无阻尼的诱导偶极模型。随后,ff12pol引入包含Thole阻尼函数的模型,从而提高了分子间相互作用能的精度
[11 ] 。
AMOEBA极化力场采用原子诱导偶极模拟极化,并使用多极矩直至四极矩表示永久静电相互作用。AMOEBA力场已经用于模拟水、离子、有机分子和蛋白体系
[44 -49 ] 。Mu等
[50 ] 采用AMOEBA力场捕捉到了
σ -hole效应。Zhang等
[51 ] 发展了适用于DNA和RNA的AMOEBA力场,并通过长达35 μs的MD模拟中进行了广泛的验证。通过模拟DNA双链、RNA双链和发夹的晶体结构,其结果与NMR确定的结构一致,均方根偏差小于2.0 Å(1 Å=0.1 nm)。
CHARMM Drude力场采用Drude振子模型描述极化以及偏离中心的电荷,从而体现电荷分布的各向异性
[6 ] 。Drude力场涵盖大多数生物分子,包括小的有机分子、蛋白/多肽、DNA
[52 ] 以及脂类
[53 ] 。Lin等
[54 ] 通过对脂肪族和芳香族卤素化分子引入异位电荷,使得卤原子表现各向异性。此外,他们还通过在Drude粒子上引入范德华参数改进Drude力场。
SIBFA
[21 ] 是一种从头算极化力场,由静电多极、短程斥力、极化、电荷转移以及色散贡献的加和构成,通过重现QM对应项来确定每一项的值。该力场首次开发是用于处理Zn
2+ 蛋白
[55 ] ,之后延伸到卤素化合物
[56 ] 以及核酸
[57 ] 中。
2.4 ABEEM极化分子力场
本课题组基于DFT和电负性均衡方法(Electronegativity equalization method,EEM)提出和发展了原子-键电负性均衡方法模型
[58 -59 ] ,随后提出了ABEEM极化分子力场(ABEEM PFF)。ABEEM PFF的特色在于将分子的单电子密度划分为原子、
σ -键、
π -键和孤对电子位点,多个位点模型体现了原子周围电荷分布的各向异性。ABEEM PFF通过求解与几何结构相关的线性方程确定目标分子的电荷分布,避免了自洽迭代过程,从而可以快速地获得电荷分布。因此,在ABEEM PFF中,原子电荷在模拟期间随着周围环境的变化而波动,有效捕捉极化效应和电荷转移现象。
2.4.1 分子和分子体系的电荷分布
分子和分子体系的电荷分布,即分布在原子或区域上的电荷是一个非常重要的物理量,对于分子内和分子间相互作用、建立分子力场以及进行分子动力学模拟都是不可或缺的。原子电荷作为描述电子密度分布的重要参数,在理论化学研究中发挥着重要作用。不同计算方法各有其物理基础和适用范围,例如波函数方法能够反映量子力学特性,而分子势能模拟方法更适用于分子力场构建。研究者们通过不断优化这些分子力场,以更准确地描述分子中的电荷分布特性,进行分子动力学模拟,从而为分子性质预测和分子设计提供可靠的理论依据。
原子电荷的计算方法根据分子区域划分方式的不同,主要分为明确边界(Sharp boundary)和模糊边界(Fuzzy boundary)两大类,具体如
图1 所示。明确边界模型基于实空间标量函数的量子化学拓扑(QCT)将分子完全划分为原子区域,将电子密度对原子区域积分可得该原子的电子数,最后将核电荷数减去电子数即得分子中的原子电荷,包括QTAIM电荷
[60 -61 ] 和PAEM电荷
[62 -63 ] 、KSpot电荷
[64 ] 等。模糊边界模型主要依据电子密度在原子间重叠,可进一步分为基于波函数,如Mulliken
[65 -70 ] 、NPA
[71 -73 ] 、Hirshfeld
[74 -80 ] 、MMPA
[81 ] 、Löwdin
[82 -85 ] 和AOIM
[86 ] 电荷等;基于电负性均衡原理计算获得的原子电荷,如EEM
[87 -88 ] 、QEq
[89 ] 和ABEEM
[59 ,90 ] 电荷等;基于拟合静电势方法获得的原子电荷,如MK
[91 -93 ] 、CHELPG
[94 ] 、ESP
[95 ] 和RESP电荷
[96 ] 等;基于实验获得的原子电荷,如根据原子极化张量获得的APT电荷
[97 -102 ] ,从红外光谱实验数据推导而得到的IR电荷
[100 ,103 ] ,基于X射线衍射获得的电荷,小分子的偶极矩也可以用来计算原子电荷
[104 ] 等;其他方法获得的原子电荷、MMFF94电荷等。
2.4.2 原子-键电负性均衡方法
ABEEM(Atom-bond electronegativity equaliza-tion method)
[23 ,105 ] 是一种基原子-键电负性均衡方法,用于计算分子中原子、键和孤对电子上的电荷分布。依据电负性均衡原理,这种方法通过求解一组线性方程来计算分子中各个原子、键和孤对电子的电荷。ABEEM特别考虑了氢键相互作用区域,在这个区域内,孤对电子和氢原子之间的相互作用使用氢键拟合函数描述。这个区域的电荷相互作用对于理解水分子之间的氢键相互作用至关重要。模型中包含了一些可调整的参数,如键电荷的位置、孤对电子的位置以及与氢原子形成氢键的孤对电子的相互作用参数。这些参数通过拟合实验数据或高精度的量子化学计算结果来确定。
图2 为胸腺嘧啶分子的ABEEM位点图,其中黄色椭圆形球体为孤对(lp)电子位点,位于双键原子共价半径上,其角度由杂化轨道类型确定。浅蓝色椭圆形球体为垂直于σ键平面的π键位点,位于双键原子共价半径上,处于双键原子的上方和下方。绿色椭圆形球体为σ键位点,位于2个键合原子之间的共价半径比处。红色、浅灰色和白色圆球分别表示氧原子、碳原子和氢原子。从
图2 中我们可以看到,胸腺嘧啶(T)共有50个电荷位点,包括原子核上的15个原子位点、15个σ键位点、16个π键位点和4个lp电子位点。ABEEM模型可以清晰地展示分子中各个位点的位置,多位点随环境变化时,电荷在所有位点之间浮动,体现了极化效应和电荷转移效应。
在ABEEM模型中,根据DFT中的有效电负性, χ i = ( ∂ E m o l / ∂ q i ) q j , … , R i j ,则每个位点的有效电负性可表示为式(6) 。
χ i = ∂ E ∂ q i R i j , q i = χ i * + 2 η i * ( q i - q i * ) + ∑ j , j ≠ i k i j q j R i j (6)
根据电负性均衡原理,分子中每个位点的有效电负性等于整个分子的有效电负性,表示为式(7) 。
χ i = χ j = ⋯ = χ ¯ m o l (7)
与此同时,所有位点的电荷之和等于分子或体系的总电荷,示于式(8) 中。
q i + q j + … = q m o l (8)
式中:
Rij 为位点
i 和
j 之间的距离;
qi 和
qj 分别是位点
i 和
j 的点电荷。当位点
i 的电荷的绝对值超过1.0e时,例如金属蛋白中的过渡金属原子以及核酸中的磷原子,需要设置价态电荷
qi * ,以使得位点
i 的相对点电荷的绝对值|
qi ′ |=|
qi -qi *|小于1.0e
[23 ] 。因此一旦确定了体系中每个位点的价态电负性
( χ i * ) 、价态硬度(
2 η i * )和静电系数(
kij ),就可以通过求解线性方程组获取不同化学环境中所有位点的电荷。
2.4.3 ABEEM极化分子力场势能函数
将ABEEM方法与分子力场结合就得到的ABEEM极化分子力场(ABEEM PFF),其势能函数表示为式(9) 。
E A B E E M P F F = E b o n d s + E a n g l e s + E t o r s i o n s + E i m p t o r s + E e l e c + E v d W (9)
式中: E b o n d s 和 E a n g l e s 分别表示为键长的伸缩振动势能和键角的弯曲振动势能,水分子的键长伸缩振动势能,使用Morse势能函数表示,其他键长伸缩振动势能则通过谐振式。 E t o r s i o n s 和 E i m p t o r s 分别表示二面角的扭曲振动势能和非共面弯曲势能项。E vdW 和E elec 分别表示范德华和静电相互作用势能,具体表达式为式(10 )~(13 )。
E b o n d s = ∑ b o n d ∈ H 2 O , c o o r d i n a t e b o n d , a n d r e a c t i o n r e l a t e d D e e - 2 α r - r e q - 2 e - α r - r e q
E b o n d s = ∑ b o n d s ∈ o t h e r k r r - r e q 2 (10)
E a n g l e s = ∑ a n g l e s k θ θ - θ e q 2 (11)
E t o r s i o n s = ∑ t o r s i o n s v 1 2 1 + c o s ϕ + v 2 2 1 - c o s 2 ϕ + v 3 2 1 + c o s 3 ϕ (12)
E i m p t o r s = ∑ i m p t o r s v 1 - c o s 2 ω (13)
式中:De 为势能最小值的势阱; α 是一个与力常数kb 有关的参数。具体表达形式为: α = ( k b / D e ) 1 / 2 ,kr 是每个键的力常数。kr 和kθ 分别表示键长伸缩振动和键角弯曲振动的力常数;r 和θ 分别表示键长和键角的真实值;r eq 和θ eq 表示平衡几何构型下的平衡键长和平衡键角;ν 1 、ν 2 和ν 3 表示二面角扭转势能项的展开力常数,ϕ 表示实际二面角值;ν 表示非共面弯曲势能项的力常数; ω 为非共面二面角真实值。
非键相互作用项,包括Lennard-Jones势能项(式(14) )和库仑项(式(15) )。
E v d W = ∑ a < b 4 f a b ε a b σ a b 12 R a b 12 - σ a b 6 R a b 6 (14)
E e l e c = ∑ i < j k i j q i q j R i j (15)
式(14) 中,原子对a 和b 的势能参数等于2个原子势能参数的几何平均,分别表示为式(16) 和(17)。
σ a b = σ a a σ b b 1 2 (16)
ε a b = ε a a ε b b 1 2 (17)
式中:a 和b 为1-2或1-3关系;系数fab 为0.0;对于1-4关系的相互作用fab 为0.5;而对于所有其他分子内和分子间相互作用fab 为1.0。
在非键势能项中除了Lennard-Jones势能项外,另一项是库仑势能项。式(15) 中,qi 和qj 表示分子中的点电荷;Rij 是点电荷qi 和qj 之间的距离;kij 是总体校正因子0.57。当i 与j 可以形成氢键时,我们采用氢键拟合函数kij = kHB ,kHB 可以表示为式(18) :
k H B ( R l p , H ) = k ' - C 1 + e x p R l p , H - r l p , H D (18)
式中:k′ 、ri,j 、C 和D 是参数,R lp,H 代表孤对电子与氢原子之间的距离。在显式水中的分子动力学模拟过程中,水分子之间存在氢键。k HB 用于处理以上水分子之间氢键的函数关系,更有效地描述水分子间静电相互作用。
2.5 分子动力学模拟
分子动力学模拟(Molecular dynamic simulation,MD),是目前应用最广泛的计算庞大复杂体系的方法。自1970年起,由于分子力学的迅速发展,人们又系统地建立了许多适用于生化分子体系、聚合物、金属和非金属材料的力场,使得计算复杂体系的结构、热力学和光谱性质的能力及准确性大为提升。分子动力学模拟是为应用这些力场及牛顿运动力学原理发展起来的计算方法。其优点为精确性高,可同时获得系统的动力学和热力学统计资料,并可广泛地适用于各种系统及各类特性的探讨。
考虑含有N 个分子的运动体系,系统的能量为系统中分子的动能与系统总势能的和。其总势能为分子中各原子位置的函数, U ( r 1 , r 2 , ⋯ , r n ) 。通常势能可分为分子间(或分子内)原子间的范德华作用(vdW)与分子内部势能(int)两部分,即表达式(19) 。
U = U v d W + U i n t (19)
范德华作用一般可将其近似为各原子对间的范德华作用的加成,分子内势能则为各类型内坐标(如键伸缩、键角弯曲…)势能的和。
系统中任一原子i 所受之力为势能的梯度见式(20) ;i 原子的加速度见式(21) ;对时间积分,可预测i 原子经过时间t 后的速度与位置,表达式为式(22) 。
F i = - ∇ i U = - i ∂ ∂ x i + j ∂ ∂ y i + k ∂ ∂ z i U (20)
a i = F i m i (21)
d 2 d 2 t r i = d d t v i = a i v i = v i 0 + a i r i = r i 0 + v i 0 t + 1 2 a i t 2 (22)
式(22) 中: r 及 v 分别为粒子的位置与速度;上标“0”为各物理量的初始值。
分子动力学计算的基本原理,先由系统中各分子的位置计算系统的势能,再计算系统中各原子所受的力及加速度,然后令式中的t = δ t ,则可得到经过 δ t 后各分子的位置及速度。重复以上的步骤,由新的位置计算系统的势能,计算各原子所受的力及加速度,预测再经过 δ t 后各分子的位置及速度。如此反复循环,可得各时间下系统中分子运动的位置、速度及加速度等资料。一般将各时间下的分子位置称为运动轨迹(Trajectory)。
3 生物分子水溶液体系的分子动力学模拟
分子动力学模拟是解析生物分子在水溶液中行为的重要工具,ABEEM-7P水模型及极化力场已成功应用到蛋白质、核酸等生物分子体系中。本章会详细介绍其应用以及ABEEM PFF在生物分子水溶液体系的分子动力学模拟方面的工作。
3.1 ABEEM-7P水模型及ABEEM极化分子力场的应用
3.1.1 ABEEM-7P水分子模型
常见的简单水分子模型
[106 -114 ] 从早期的ST2和ST4模型逐渐发展到简单点电荷力场(Simple point charge, SPC),以及可转移的分子间势能函数(Transferable intermolecular potential, TIP
n P),这些模型均利用刚体水分子结构,即固定的键长和键角。水分子间的势能面利用van der Waals相互作用,及三点、四点或五点间的静电相互作用来描述。如
图3 中常使用的三点、四点和五点的简单水分子力场模型。
这些模型均利用刚体水分子结构,同时以固定电荷为参数,无法描述原子(或分子)周围电场改变而产生的静电相互作用的变化。而我们建立了一种新的计算水分子体系的模型,可转移的分子间七点浮动电荷非刚体模型(a Transferable Intermolecular Potential Seven Points Approach Including Fluctuating Charge and Flexible Model, ABEEM-7P)
[115 ] 。ABEEM-7P模型利用更多的位点来计算水分子间的静电相互作用,
图4 为ABEEM-7P水分子模型的单体水结构。一个单体水分子划分为7个位点,3个原子核严格定位于氢、氧原子的几何坐标。氧原子通过
sp 3 杂化形成四面体电子分布,2个σ键(O—H键)电荷密度极大值位点沿σ键分布,其精确位置由Pauling共价半径比例(氧原子0.66 Å,氢原子0.28 Å)决定。2个σ键间形成104.5°的键角,这一角度源于氧原子孤对电子对成键轨道的空间排斥效应。此外,两对孤对电子(lp)位点以109.47°的立体角对称分布在氧原子核的周围,距离氧原子中心0.74 Å(约为氧原子范德华半径的40%),其空间构型严格遵循
lp 间排斥能最小化原则。
在ABEEM-7P模型中,键长的伸缩振动使用的Morse势能函数。由于液态水体系中氢键的形成使得H—O键有不同的存在形式,如OH- 和H3 O+ ,Morse势能函数能够从平衡键长到键的解离更大范围描述键的伸缩振动,对详细分析水分子体系的振动模型极为重要。
ABEEM-7P模型中O—H键的平衡键长
r eq 和H—O—H的平衡键角
θ eq 分别被设置在实验值0.9572 Å和104.52°。ABEEM-7P模型的力场参数列于
表1 中,其中氧-氢键的离解能
D 为529.6 kcal/mol,角力常数
kθ 为34.05 kcal/mol·(°)
2 。α是和键长的力常数相关的参数(
α = ( k b / 2 D ) 1 2 ,
kb 是键的力常数),同时给出了TIP4P-FQ
[116 ] 、TIP5P
[107 ] 和POL5
[117 ] 模型的力场参数。
使用ABEEM-7P水分子模型计算得到的单体水电荷分布和偶极矩。相比于其他水分子力场,ABEEM-7P水分子模型充分考虑了构象的变化,给出了所有原子、键和孤对电子的电荷分布。同时ABEEM-7P模型得到的偶极矩1.8551与实验值1.855
[115 ] 符合得非常好,这与各位点合理的部分电荷密切相关,因此精确的描述静电相互作用可以准确地描述水分子的偶极矩。
3.1.2 水团簇
水分子团簇是众多理论和实验研究的热门课题之一,研究水分子团簇的一个重要目的是希望对水溶液行为做出分子水平的定量描述。高水平的水分子团簇研究可望为深层次揭示化学、生命科学和信息科学等领域的本质问题提供有力工具。水分子团簇结构与功能研究的重要理论意义和广泛的应用前景已使其成为活跃的科学前沿研究热点之一。
从结构上看,水分子团簇可以分为两种类型,一种是二维结构(如环状结构),另一种是三维结构(如笼状、棱状或者立方状结构)。需要进一步说明的是,在水分子团簇中有4种类型的单体,即1个质子给予体和1个受体(“da”型),1个给予体和2个受体(“daa”型),2个给予体和1个受体(“dda”型),2个给予体和2个受体(“ddaa”型)。
图5 给出了(H
2 O)
n (
n =2~10)的各种低能结构。在七聚水的5种低能结构中,缺少一个角的类似于立方状的结构7A是最低能量结构,能量次低结构是7B,这和从头计算
[119 -120 ] 所报道的结果是一致的,实验上的光谱数据
[121 ] 也证实了7A和7B这两种结构的存在。通过观察,7A和7B可以从具有对称性S4的八聚水中移走一个“dda”或者“daa”型的单体水分子而得到;同时也可以通过在六聚水棱状(prism)结构之上或之下的三角形基准面上插入一个“da”型的水分子而得到。值得注意的是,
图5 中的所有结构(除了结构D2d)通过反映之后都能够得到它们各自的等能异构体
[120 ] ,这些异构体和先前的结构是不能重叠的,即互为手性对称。7A有10个氢键,结合能是-230.25 kJ/mol。7B也有10个氢键,其结合能是-227.65 kJ/mol,比7A稍微高出了2.60 kJ/mol,与Kim等
[119 ] 用MP2/6-311++G**//HF/6-311++G**方法、且考虑了50%的基组叠加误差校正(BSSE),以及用MP2/6-311++G**//B3LYP/6-311++G**方法、也加入50%的BSSE所得到的两者结合能之差约为1.90 kJ/mol的结果相差甚微,也相似于Sathyamurthy等
[122 ] 在HF/6-31G(d,p)//HF/6-31G(d,p)水平下得到的结果。双环状结构R4和R5都可以被看作是在六聚水书状(book)结构的基础上插入一个水分子而得到的,它们有相似的结合能和相同的氢键个数。
该模型预测了七聚水中各种低能结构的偶极矩和结构参数。对于偶极矩,该模型的计算结果与Kim等
[119 ] 在MP2/TZ2P++//B3LYP/6-311++G**水平下、以及我们在HF/6-31G(d,p)//HF/6-31G(d,p)水平下的计算结果相比,都呈现了很好的一致性。除此之外,ABEEM-7P模型已经应用于更大的水分子团簇(H
2 O)
n (
n =20, 24, 28, 32, 35)的研究中,结果表明,对于更大的水分子团簇,在团簇的几何构型和结合能等方面,ABEEM极化分子力场模型也能得到与从头计算或其他半经验方法相一致的结果。以上结果将有助于加深对中等大小的、较大的乃至更大的水分子团簇的认识,并对水溶液性质的了解有重要的参考价值,也为进一步验证ABEEM极化分子力场模型的合理性、参数的正确性和可转移性提供了可靠的依据。
3.1.3 多肽和蛋白质
生物分子的构象分布研究,对了解生物分子的功能是非常重要的。稳定构象的正确搜索是挑战性的任务。2015年,Xu等
[123 ] 使用ABEEM PFF用于搜索Val-dipeptide和Val-tripeptide的稳定构象,分别得到6种和34种稳定构象,这与从头算方法B3LYP/6-311++G(d,p)的结果是一致的,而OPLS/AA-FCFF和AMBER-FCFF 99sb力场只能分别搜索到较少或更少的稳定构象。
2011年,Zhao等
[105 ] 将ABEEM模型与GB方法结合计算的18种氨基酸侧链(除GLY, PRO)小分子的溶剂化自由能和其他方法计算结果
[124 ] 与实验结果的偏差进行了统计,我们的溶剂化自由能相对与实验结果偏差是最小的,为0.17 kcal/mol,说明ABEEM方法十分适合GB模型用于计算溶剂化自由能。
为了进一步测试ABEEM模型中用于计算连续溶剂中电荷的校正系数,我们选择20种二肽共有120种稳定构象,在PCM连续溶剂模型下,用MP2/AUG-CC-PVDZ和B3LYP/AUG-CC-PVDZ量子力学方法计算了这120个构象的溶剂化自由能,分别分析了ABEEM/GB计算二肽的溶剂化自由能的线性关系。其中,ABEEM/GB与MP2/AUG-CC-PVDZ方法的线性相关系数为0.9951,线性斜率为1.0625;ABEEM/GB与B3LYP/AUG-CC-PVDZ方法的线性系数为0.9945,斜率为0.9783。很明显,我们的方法预测二肽溶剂化自由能的能力与从头算方法和DFT方法都是非常吻合的,这为之后的生物系统模拟提供了有力的支持。
2011年,Chen等
[125 ] 将ABEEM极化分子力场与LIE方法结合得到了ABEEM/MM LIE方法,用于估算蛋白质-配体绝对结合自由能。通过使用24个HIV-1蛋白酶-抑制剂复合物的训练集对该方法的参数进行校准,得到了双参数模型和三参数模型24个受体-配体复合物的相互作用能和溶剂化能。结果显示,三参数模型比双参数模型具有更高的预测精度。其中,RMSD和平均绝对偏差分别为0.77、0.70 kcal/mol和0.55、0.56 kcal/mol。为了避免数据过度拟合,我们在模型上进行了留一交叉验证,交叉验证结果显示三参数模型的
q L O O 2 和
S P R E S S L O O 分别为0.82和0.78 kcal/mol,双参数模型和三参数模型的多重相关系数
r 2 分别为0.67和0.87。同时第三个参数
Δ G t r , r o t =8.67 kcal/mol在实验得到的7~11 kcal/mol的范围内
[126 -127 ] 。双参数模型和三参数模型得到24个复合物的自由能与实验数据的线性相关性分别为0.91和0.93,进一步证明了三参数模型预测精度更高。
为了进一步验证模型的预测能力,将模型用于7个不同结构的HIV-1蛋白酶-抑制剂复合物,得到的结果也是令人满意的,三参数模型计算的结合能结果与实验值的RMSD为1.24 kcal/mol。目前此方法相对结合自由能预测的普适性较好,但参数需要根据不同的靶标重新校准,可通过分析回转半径和氢键距离作为结合能评估的辅助指标。
2014年,Yang等
[23 ] 建立了金属蛋白中所有过渡金属原子的ABEEM参数。通过引入参考电荷用来处理过渡金属可能存在的较大电荷;通过线性回归和最小二乘法校准了价态电负性(Valence state electronegativity)和价态硬度(Valence state hardness)参数。选择700个含有过渡金属原子的生物分子为模型分子,优化了电荷相关参数,并与HF/STO-3G水平的从头算结果进行对比验证。以Mo不同价态下的25个金属蛋白为例,得到模型分子的ABEEM电荷与从头算方法(HF/STO-3G)计算电荷的线性关系,所有的斜率都接近1.0,相关系数
R 都超过了0.97,说明了2种电荷具有很好的线性关系。以5种金属蛋白为例,它们的PDB编号分别为3bdja、1ti4、1vlfj、1g8ja和1g8ka,对比了ABEEM极化分子力场优化后的结构与从头算MP2/6-311G方法优化的结构,结果2种方法的结构非常接近,其中键长的RMSDs分别为0.0396、0.2158、0.2008、0.1309和0.0238 Å,键角的RMSDs分别为7.93°、8.91°、9.25°、6.96°和6.06°,说明ABEEM极化分子力场可以用于优化金属蛋白的几何结构,这为后续金属蛋白的动力学模拟提供了基础。
3.1.4 核酸
DMP是核酸分子中磷酸二酯键重要的模型分子。2009年,Wang等
[128 ] 使用ABEEM方法构建了含有35个位点的DMP模型,研究了气相DMP以及DMP-(H₂O)
ₙ (
n =1~3)团簇的性质。与QM方法相比,ABEEM方法优化DMP结构的键长、键角和二面角的RMSDs分别为0.027 Å,2.4°和1.1°。对DMP-(H₂O)
ₙ (
n =1~3)团簇的氢键结构与相互作用能进行分析。其中,ABEEM和QM方法优化的DMP-H
2 O二聚体的氢键角和氢键长的RMSDs分别为1.2°和0.020 Å;使用ABEEM方法得到了11种DMP-(H₂O)
ₙ (
n =1~3)团簇的相互作用能与MP2/6-311++G**方法的计算偏差≤1.94 kcal/mol,说明ABEEM PFF能够准确处理DMP-H
2 O复合物的氢键相互作用。
通过对ABEEM PFF势函数的参数化,对显式水分子环境下的DMP阴离子水溶液进行了分子动力学(MD)模拟,发现DMP阴离子的非酯化氧原子周围形成了明确的第一水合层,其平均配位数为5.2
[128 ] 。特别值得注意的是,我们发现一个水分子的2个氢原子可同时与DMP阴离子的2个非酯化氧原子形成双氢键,这是ABEEM PFF与AMBER ff99
[13 ] 力场在MD模拟中最显著的差异,同时也证明浮动电荷力场比传统固定电荷力场能提供更丰富的信息。
核酸(如DNA)的生物功能与其构象密切相关。2009年,Wang等
[129 ] 得到了ABEEM PFF糖-磷酸骨架和糖苷键的二面角扭转势能参数,并通过氢键拟合函数增强对氢键相互作用的描述。ABEEM PFF能够很好描述核酸碱基分子的结构、相互作用能和偶极矩。通过ABEEM方法优化的甲基化碱基(腺嘌呤A、鸟嘌呤G、胞嘧啶C、胸腺嘧啶T)非氢原子键长与晶体学数据的RMSD为0.005~0.006 Å,与CHARMM27
[20 ] 力场结果相当。氢原子相关键角的RMSD为0.4°~3.2°,其中鸟嘌呤因氨基平面化假设偏差较大(3.2°),其他碱基偏差较小(如胸腺嘧啶仅0.4°),说明ABEEM PFF能够很好描述核酸碱基分子的结构。使用ABEEM PFF计算的甲基化碱基偶极矩与实验值及高精度理论方法
[130 -132 ] 吻合良好,相对于MP2理论的RMSD为0.27
D ,优于CHARMM27
[61 ] 力场的0.84
D 。同时,ABEEM PFF能准确复现核酸体系中关键二面角的势能面
[129 ] ,例如二面角
γ 在C3'-endo的3个极小值位置分别在43.1°、191.0°和298.4°,这与从头算方法得到的44.1°、191.1°和298.3°基本一致。通过对糖苷键
χ 的势能面分析,显示反式构象(anti)更加稳定,这与晶体结构是一致的。
通过对Uracil-(H
2 O)
n (
n = 1~7)团簇的38个结构及其能量的研究
[31 ] ,结合高精度的量子化学计算(MP2/6-311++G**),探讨了水分子与尿嘧啶间的氢键相互作用。通过对比ABEEM PFF和从头算方法的相互作用能,由
图6 可知,得到的所有能量的RMSD为0.96 kcal/mol,线性相关系数为0.997,说明ABEEM PFF能够很好的模拟碱基与水的相互作用。
此外,2014年,Liu等
[133 ] 以B3LYP/6-31+G*水平的量子力学计算结果为基准,研究了8-氧代腺嘌呤(8-oxo-A)碱基对的电荷分布、能量、结构。与量子力学结果相比,ABEEM PFF计算的电荷分布、偶极矩和相互作用能的线性相关系数分别达到0.99、0.96和0.96,键长与键角的均方根偏差(RMSDs)分别仅为0.005 Å和0.5°。在此基础上,2021年又研究了修复酶hOGG1通过氢键和π-π堆叠作用识别氧化碱基(8-oxo-G)的分子机制
[30 ] ,重点揭示了水分子在识别过程中的动态介导作用。hOGG1中的Gln315通过2个水分子(A和B)的连续交换与8-oxo-G形成氢键,Phe319和Cys253从底部π-π堆叠在平面的两侧,形成三明治结构。电荷极化效应提供了一个重要的信号来驱动水分子的交换,并保持酶对氧化碱的识别。这2个介导水分子一起通过氢键相互作用拉动Gln315识别8-氧代-G。溶剂水分子C的电荷极化效应(电荷绝对值增加约0.2 e)驱动A/B的交换,维持酶对氧化碱基的识别。通过对比量子力学(QM)计算结果与不同力场的模拟结果,验证了ABEEM PFF在几何结构、电荷分布和氢键能量方面的准确性。对比ABEEM PFF、OPLS/AA和AMBER OL15力场
[134 ] 与QM(M06-2X/6-311++G**方法)优化模型分子的几何参数(键长、键角、二面角)的RMSDs,ABEEM PFF的RMSD最小(键长0.017 Å,键角2.20°,二面角2.26°),表明其几何结构拟合最优,AMBER OL15在二面角模拟中偏差较大(RMSD达9.50°),导致甲基构象和平面结构失真。ABEEM PFF的偶极矩与QM结果的平均绝对偏差(AAD)为1.45
D ,显著优于OPLS/AA和AMBER OL15的3.70
D 和2.95
D (偏差为2~3倍)
[134 ] 。ABEEM PFF计算的氢键能量与MP2/aug-cc-pVDZ方法的AAD为2.61 kcal/mol,RMSD为3.30 kcal/mol,明显优于OPLS/AA和AMBER OL15的AAD和RMSD
[134 ] ,其中,OPLS/AA的AAD为10.35 kcal/mol,RMSD为15.65 kcal/mol,AMBER OL15的AAD为13.50 kcal/mol,RMSD为25.52 kcal/mol。说明ABEEM PFF能准确描述分子的几何结构、氢键作用区域的静电相互作用以及电荷分布。
3.2 ABEEM PFF在生物分子水溶液体系的分子动力学模拟实例
3.2.1 GA 88和GB 88的分子动力学模拟
2013年,Wang等
[135 ] 选用G
A 88 (PDB编码: 2JWS)和G
B 88(PDB编码:2JWU)2个分子,使用ABEEM PFF和OPLS-AA固定电荷力场方法进行了温度为295 K的水溶液和真空分子动力学模拟。比较2种方法得到的两个蛋白质的结构与实验结构的均方根偏差,分析了2种方法得到回旋半径、氢键及电荷分布情况。结果表明,ABEEM PFF得到的GA88和GB88中非氢原子的RMSDs均小于OPLS-AA固定电荷力场;回旋半径模拟很好地体现了蛋白质的“电致紧缩”现象;ABEEM PFF模拟的蛋白质GA88和GB88的氢键平均键长与NMR结构中相应的氢键平均键长的RMSDs分别为0.012和0.008 nm,OPLS-AA力场的结果则分别为0.012和0.011 nm,两者基本一致。ABEEM模型不但考虑了原子区域的电荷,而且加入了化学键和孤对电子区域电荷,各区域电荷能够随体系环境的改变而及时计算变化后体系中各位点的电荷。而利用OPLS-AA固定点电荷力场进行模拟时,无论什么环境下,所有原子的电荷都不变。因此,ABEEM PFF能够更好地模拟蛋白质的结构并体现蛋白质和周围水分子的极化效应。
3.2.2 硼酸-β -内酰胺酶的分子动力学模拟
硼酸作为
β -内酰胺酶抑制剂,在治疗生物感染和肿瘤中具有潜力。2020年,Lu等
[136 ] 构建了适用于硼体系的ABEEM PFF,弥补了现有力场缺少硼参数的空白。ABEEM PFF的优势在于可以提供合理的电荷分布、电荷转移以及对氢键相互作用的特殊处理。选取5类共47硼酸-
β -内酰胺酶为模型分子确定并验证了ABEEM PFF的相关参数。通过与量子力学(QM)计算结果对比,模型分子2种方法计算电荷的线性相关系数超过0.97;两种方法计算的偶极矩AAD为0.97
D ,RMSD仅为1.13
D ,说明ABEEM PFF计算的电荷分布是合理的。将2种方法优化的几何结构进行对比,确定了ABEEM PFF中Morse参数、键角参数和二面角参数的合理性。研究了10种硼酸
β -内酰胺酶模型分子和10种硼酸-
β -内酰酶与水复合物的结合能,ABEEM PFF的计算结果与MP2方法的结果具有良好的相关性。在此基础上,对5种大型硼酸-
β -内酰胺酶复合物(BZB-1C3BJ、CTB-1FSWJ、NPB-1KDSJ、SM3-1MY8J和0NB-4E3IJ)进行了MD模拟,模拟时间为4.0 ns。我们对MD模拟的结构与PDB的初始结构进行对比,其中,5种结构键长的RMSDs分别为0.041、0.039、0.041、0.030和0.044 Å;键角的RMSDs分别为3.89°、3.03°、3.99°、3.73°和3.66°。上述结果说明ABEEM PFF可以模拟得到硼酸-
β -内酰胺酶的稳定结构。由于二面角扭转的势垒较低,因此很容易受到周围环境的影响而导致与初始结构发生较大变化,这也反映了ABEEM PFF能够很好地处理实际环境中产生的极化效应。
3.2.3 镁离子蛋白和镁离子水溶液的分子动力学模拟
镁离子是重要的生物金属离子,会参与许多重要的生理过程。2025年,Zhang
[137 ] 等建立了镁离子的ABEEM PFF的成键模型(ABEEM-BM)和非键模型(ABEEM-NBM),用于模拟镁离子-蛋白质复合物和水溶液中的表现。截取40个PDB中的镁蛋白片段,保留镁离子、配体氨基酸侧链及水分子,通过限制性优化获得稳定结构。使用ABEEM-BM、ABEEM-NBM、OPLS-AA
[38 ] 、AMBER99
[138 ] 和CHARMM22
[35 ] 力场,对水溶液中的40个镁蛋白片段进行了10 ns的MD模拟。与PDB结构相比,通过MD模拟的RMSDs(键长/角度),ABEEM-BM,ABEEM-NBM,OPLS-AA,AMBER99和CHARMM22的键长的平均RMSDs为0.088、2.490、2.697、2.831和2.594 Å。平均角度的RMSDs分别为5.99°、24.38°、24.96°、25.86°和24.01°。ABEEM-BM的键长和键角RMSD明显小于ABEEM-NBM,OPLS-AA,AMBER99和CHARMM22力场,说明ABEEM-BM能够得到镁蛋白分子的稳定构象。
使用ABEEM-BM和ABEEM-NBM进行镁离子水溶液的MD模拟,研究了镁离子在水溶液中的动态性质。对径向分布函数(RDF)的分析结果如下,ABEEM-BM的Mg—O第一峰位置出现在2.050 Å,与实验值(2.090 Å)
[139 ] 和AIMD(2.100 Å)
[140 ] 结果一致,而ABEEM-NBM出现在1.850 Å。2种方法的获得的Mg—O第一峰积分都是6,说明水溶液中镁离子的配位数为6,与实验值和AIMD结果相匹配。ABEEM-BM和ABEEM-NBM 2种方法水交换速率常数分别为4.7×10
5 和5.0×10
6 s
-1 ,ABEEM-BM更接近实验值(5.3×10
5 ~6.7×10
5 s⁻¹),相对于实验值误差小于microMg(TIP3P)
[141 ] 。镁离子第一个水合层内的水分子的平均停留时间(
τ )实验值为1.0×10
-6 s
[142 ] 。与ABEEM-NBM(2.0×10
-7 s)相比,通过ABEEM-BM(2.1×10
-6 s)计算得出的
τ 更接近实验数据。这些结果表明,ABEEM-BM在维持计算效率的同时,提高了镁离子水溶液动力学模拟的准确性。
图7 描述了ABEEM模型动态捕捉镁离子与配体间的电荷转移过程,镁离子部分电荷区间为1.07~1.24|e|(ABEEM-BM),与QM结果(1.06~1.42|e|)高度吻合。这项工作可为后续利用ABEEM PFF研究水合镁蛋白的结构与性质奠定了基础。
3.2.4 血红素蛋白酶催化苯乙烯氨羟化水解生成对映异构体
在设计的血红素蛋白酶调控催化反应的活性和选择性中,空间结构效应和静电效应扮演着关键角色,而极化力场能有效地描述这种调控机制中由结构变化引起的静电效应的改变
[143 ] 。
2024年,Huang等
[144 ] 聚焦于设计的细胞色素
c (P450)在催化碳自由基Fe卟啉复合物水解过程中的重要功能。利用模型分子拟合氨羟化反应相关的ABEEM PFF的电荷参数,通过ModSeminario
[145 ] 和扫描方法,确定了键长伸缩和键角弯曲参数,并将Fe与轴向配体间的相互作用引入Morse势能函数描述。通过构象搜索,确定了活性中间体的结构。诱变的残基通过vdW相互作用和100S-101P残基的氢键相互作用对活性中间体施加了限制,影响其结构的朝向。同时,活性中间体与水分子之间的氢键相互作用共同促进了
R 型产物的形成。此外,QM/MM-ABEEM PFF和FEP计算表明,在3Rts结构中,反应复合物与相邻氨基酸残基之间能形成稳定的三明治状结构,相比于3Sts中的拥挤排列,3Rts显示出更高的稳定性(
图8 )。
R 型对映体的自由能势垒比
S 型对映体低约10.0 kcal/mol,这使得
R 型对映体的形成更为有利
[144 ] 。这些研究为设计细胞色素
c 在苯乙烯氨羟化反应中具有优异的立体选择性提供了机制上的见解。
3.2.5 DNA的变异和修复机制
生物体内的修复酶对受损DNA具有修复作用,其本质是酶与受损DNA发生的酶催化反应,揭示该反应的微观机制将对相关疾病的治疗提供帮助。2020年,Liu等
[146 ] 通过ABEEM PFF对4种不同修复酶切除受损碱基的20多条反应路径进行了水溶液体系的分子动力学模拟。所有反应均为糖苷键断裂协同质子转移的协同SN
2 反应。ABEEM PFF定义了反应价态电负性函数,精准地模拟了反应过程中,参与化学键形成和断裂原子及过渡态附近相关原子电荷的强烈变化;通过变化的局域电荷守恒条件模拟出由活性残基到受损碱基的约-1.0e的电荷转移过程。分子动力学模拟结果表明,质子化碱基是更好的离去基团。酶切除质子化碱基反应的活化能比切除中性碱基反应的活化能普遍低10 kcal/mol。对于受损鸟嘌呤
[30 ] 和受损腺嘌呤
[147 ] ,修复酶采用双功能机制进攻受损碱基切断糖苷键;对于受损胞嘧啶
[148 ] 和受损胸腺嘧啶
[149 ] ,修复酶采用单功能机制进攻受损碱基切断糖苷键,
图9 展示了这一过程。ABEEM极化反应力场为探究水溶液环境中DNA糖基化酶切除损伤碱基反应提供了新的理论方法。我们将对整个蛋白和DNA体系的水溶液进行模拟,探讨远端残基对酶识别的影响,揭示酶特异性识别的本质。
3.2.6 鸟嘌呤自由基阳离子(G·+)去质子化机制
G·+的快速去质子化会阻碍DNA链中空穴(正电荷)的远距离传输,研究其机制有助于设计高导电性DNA结构
[150 ] 。2024年,Wang等
[150 ] 在3H
2 O和9H
2 O模型中,通过QM和QM/MM(ABEEM)方法详细讨论了电荷分布、电荷转移、几何结构和活化能(
图10 )。
由
图10 可知,QM方法获得的3H
2 O和9H
2 O模型的质子转移距离和活化能分别为0.469 Å、36.6 kJ/mol和0.391 、26.0 kJ/mol。3H
2 O模型中G·+的相对能量远大于9H
2 O模型,这表明后者的质子转移更容易,同时,9H
2 O模型中第二水合层的存在降低了势能势垒,促进了质子转移。使用QM/MM(ABEEM)方法进一步模拟了9H
2 O模型中G·+的微观脱质子化机制。利用价态电负性(χ*)分段函数和引入电荷局部守恒条件,该方法可以精确模拟极化和电荷转移的影响。我们拟合了QM和QM/MM(ABEEM)方法计算的电荷分布之间的线性关系,相关系数超过0.97,说明2种方法之间的电荷转移是一致的。MD模拟表明,脱质子化活化能为(16.3±0.8) kJ/mol,与实验结果((15.1±1.5) kJ/mol)
[151 ] 一致。上述结果表明QM/MM(ABEEM)方法在模拟质子转移中的成功应用。G·+在游离碱中的脱质子化速率(1.5×10
7 s
-1 )快于双链DNA(10
6 ~10
7 s
-1 )。此外,质子的远距离转移对应于高能势垒。在设计DNA电子器件时,建议避免加入自由碱基,优先考虑具有G·+远距离质子转移的分子。上述机理对电子器件结构的合理设计具有指导意义。
3.2.7 脂质双分子层的分子动力学模拟
细胞膜的结构以脂质双层为基本框架,对多种生物过程至关重要。为深入理解其特性及相互作用机制,学界已对其开展了广泛研究。2025年,Wang等
[152 ] 构建了用于PC脂质的ABEEM PFF(
图11 ),合理地再现了脂质双分子层的关键实验特性,包括体积、厚度、面积、核磁共振(NMR)有序性脂质双层的参数和电子密度。
对于体积和面积,DMPC模拟值分别为60.5 Å
2 和1230 Å
3 ,这与实验值
[153 ] 60.6 Å
2 和1229 Å
3 是一致的;相对于实验结果,使用ABEEM PFF优化后的碳原子氘序参数明显优于Drude PFF
[154 -155 ] ;双层上的电子密度分布(EDPs)与实验X射线散射数据表现出很好的相关性,RMSD为0.05 e/Å
3 。在主体水区域预测的偶极电势为0.403 V,这是反映脂质极化率的一个重要特征。ABEEM PFF有效地捕捉了电荷波动,并提供了更准确的电荷分布表示。根据这一特征,我们使用无监督机器学习方法研究了脂质中P原子的电荷分布与其构型之间的相关性,揭示了P-N矢量相对于膜法线方向的变化具有电荷依赖性。此外,我们的分析表明,脂质双层中P原子的电荷会影响水分子的动态行为,P电荷越高,水分子的聚集程度就越高。推进这一模型需要开发算法、软件和应用方面的努力。当前,非并行化的ABEEM PFF模拟的计算成本确实高于经典力场(FF)的成本。但使用ABEEM-DBSS
[146 ] 方法对于所有考虑的脂质双层(14,112~16,704个原子,包括72个脂质和2160个水分子),ABEEM PFF代码在AMD Ryzen 9 9950X上获取1 ns数据的时间大约是AMOBEA 09的1/4。这说明程序并行化对于将ABEEM PFF应用于复杂脂质膜和膜蛋白的模拟是十分必要的。
4 结论与展望
生物分子水溶液的研究在生命科学、医学和工业等领域具有重要意义。分子动力学模拟是研究生物分子在液态环境中结构与动态行为的重要工具。本文以ABEEM PFF为主要研究对象,详细介绍了ABEEM-7P水模型及ABEEM PFF在生物分子体系中的应用,通过一些具体实例介绍了ABEEM PFF在生物分子水溶液的分子动力学模拟中的应用。
(1)ABEEM模型含有原子、键以及孤对电子位点,多个位点模型体现了原子周围电荷分布的各向异性,特别考虑了氢键相互作用区域。通过求解与几何结构相关的线性方程确定目标分子的电荷分布,避免了自洽迭代过程,从而可以快速地获得电荷分布。ABEEM-7P水模型对理解水分子之间的氢键相互作用至关重要,因此能够精确地描述水溶液的静电相互作用。
(2)ABEEM PFF不仅可以很好地描述生物分子在水溶液中的结构和变化,同时还可以计算其溶剂化能和自由能等性质。ABEEM PFF在模拟期间原子电荷随着周围环境的变化而波动,能够很好地描述分子体系的极化效应和电荷转移。目前已经成功应用到水团簇、蛋白质、核酸和脂质等生物分子的水溶液体系。
综上所述,本文通过对大量的相关研究工作进行总结,证明ABEEM-7P水模型和ABEEM PFF可以用于生物分子水溶液的分子动力学模拟,并得到了比较好的结果。目前,仍存在一些急需解决的问题,算法和软件的开发、新模型以及应用研究对我们模型的发展是至关重要的。在未来的工作中,将机器学习方法与ABEEM模型结合也将是重要的发展方向。