基于化學鍵偶極的精氨酸分子力場參數的建立

劉 琦, 祝佳怡, 姜笑楠, 孫長亮, 郝 強, 王長生

(1.遼寧師范大學 化學化工學院,遼寧 大連 116029; 2.沈陽化工大學 分析測試中心,遼寧 沈陽 110142)

在細胞代謝、信號傳導以及基因表達等關鍵生命過程中,精氨酸顯現出不可或缺的作用. 作為蛋白質的構建元素之一,精氨酸參與了蛋白質的構建,氨基酸及其衍生物的生物合成過程[1],同時也參與了尿素循環的調控[2]. 精氨酸在細胞質和細胞核內的生物合成過程中也具有重要的功能. 此外,精氨酸在胍基乙酸以及隨后形成的肌酸合成過程中充當了重要的供體. 這種多重作用使得精氨酸成為維持蛋白質結構穩定性和功能活性的關鍵成分. 精氨酸在染色體包裝、DNA修復以及其他細胞過程中的參與[3],使其備受生物學研究關注.

然而,精氨酸分子的復雜性遠不止于其在生物學中的角色. 分子內和分子間的相互作用很大程度上受到其胍基上電荷分布和立體構型的影響. 這些電荷分布在精氨酸分子與其他分子之間的氫鍵、離子鍵以及范德華力等相互作用中扮演著至關重要的角色. 精確描述這些相互作用對于理解精氨酸在蛋白質結構的穩定性、酶活性以及其他生物學功能方面的作用至關重要[4].

分子力場方法作為分子模擬和計算化學領域的基石,為揭示生物分子的結構和相互作用提供了關鍵支持. 然而,現有的力場方法在涉及帶電體系的精確描述方面仍存在限制. 特別是在描述如胍基上電荷分布這類復雜氨基酸的時候,現行的力場參數無法提供令人滿意的結果. 因此,開發出適用于精氨酸分子的精準力場參數成為當前亟待解決的任務. 通過建立一個基于化學鍵偶極的精氨酸分子力場參數,將能夠更加準確地模擬精氨酸胍基上電荷分布以及與其他分子間的相互作用,從而有助于更加深入地研究和理解生物分子中精氨酸的功能.

1 化學鍵偶極的分子力場方法

為了精確預測精氨酸分子,尤其是帶電氨基酸側鏈部分與水分子之間的相互作用,本文對質子化的精氨酸側鏈與水分子組成的團簇體在NPT系綜下進行2 ns的分子動力學模擬. 以軌跡中精氨酸側鏈分子及其周圍的兩個水分子組成的復合物為基礎,在B3LYP/6-31+G(d,p)水平下進行了結構優化. 以一個質子化的精氨酸側鏈分子與兩個水分子形成的復合物為研究對象,構建了一個包含15個穩定三體復合物(S1～S15)的訓練集(如圖1所示). 此外,還在考慮基組重疊誤差校正的情況下,采用MP2方法結合aug-cc-pVTZ基組,獲得了訓練集中各構象的二體和三體相互作用能. 二體和三體相互作用能分別由式(1)、式(2)和式(3)、式(4)計算得到.

圖1 15個由精氨酸側鏈與水形成的三聚體復合物Fig.1 15 Trimeric complexes formed by 15 arginine side chain interacting with water molecules

V2b(i,j)=E(i,j)-E(i)-E(j),

(1)

(2)

V3b(i,j,k)=E(i,j,k)-E(i,j)-E(i,k)-E(j,k)+E(i)+E(j)+E(k),

(3)

(4)

基于課題組前期的研究基礎,本文采用了靜電作用函數、范德華作用函數和軌道作用函數構建了一種適用于質子化精氨酸分子的分子間相互作用勢函數. 與傳統的力場方法不同,這個勢函數的靜電相互作用部分是基于可極化的化學鍵偶極矩建立的,其具體函數形式如式(5)所示:

(5)

(6)

其中,δμA-B代表了分子中A-B鍵的誘導偶極矩,μ0,A-B是無外部作用時A-B鍵的固有偶極矩,qind,A和qind,B分別是采用AM1方法計算得到的復合物的A-B鍵上A原子和B原子的電荷,qind,A,0和qind,B,0分別是無外部作用下分子中A-B鍵上A原子和B原子的電荷.

在前期研究工作的基礎上[5-8],本文對精氨酸側鏈相關參數,尤其是胍基部分進行了參數化. 表1給出了本文確定的相關參數.

表1 本文所確定的參數Table 1 Parameters determined in this work

基于上述參數,計算了訓練集中的15種構象的二體和三體作用能,如表2所示. 表中列出了采用MP2/aug-cc-pVTZ和本文方法計算得到的二體作用能和三體作用能的數值,以及本文方法與MP2/aug-cc-pVTZ的絕對誤差. 該結果顯示,重新擬合的參數能夠準確計算這些復合物的二體和三體作用能.

表2 采用不同方法計算15個質子化的精氨酸側鏈與水分子復合物的二體作用能、三體作用能及相應絕對誤差Table 2 Two-body interaction energies,three-body interaction energies,and corresponding absolute errors calculated using different methods for 15 protonated arginine side chain-water complexes kJ·mol-1

2 應用

利用前述優化后的偶極-偶極作用模型式(5)和式(6)和表1中提供的參數,對一系列含有質子化精氨酸側鏈的水分子團簇[Arg(+) side chain·(H2O)n,(n= 8,10,12,15,18,20)]的二體和三體作用能進行了研究,并與高精度的MP2/aug-cc-pVTZ水平結果進行比較. 一系列含有質子化精氨酸側鏈的水分子團簇的結構均取自分子動力學模擬. 分別以軌跡中精氨酸側鏈分子及其周圍的8、10、12、15、18和20個水分子組成的復合物為基礎,在B3LYP/6-31+G(d,p)水平上進行結構優化,得到了如圖2所示的6個質子化精氨酸側鏈與水分子形成的穩定團簇結構.

圖2 質子化的精氨酸側鏈分子分別與8、10、12、15、18和20個水分子形成的團簇Fig.2 Clusters formed by protonated arginine side chain with 8,10,12,15,18,and 20 water molecules

從測試團簇分子結構上講,6種不同尺寸的水分子團簇中含有大量的分子間氫鍵作用. 不同于中性團簇分子,質子化的精氨酸側鏈分子與水分子間的氫鍵作用更強,構建的穩定團簇中,N-H…O型氫鍵的平均鍵長僅為0.211 nm. 其中,最短氫鍵鍵長為0.184 nm,最長氫鍵鍵長為0.293 nm.

二體和三體相互作用在分子間相互作用研究中扮演著至關重要的角色,其重要性不言而喻. 事實上,這兩種相互作用的貢獻在分子之間的總相互作用中可占95%以上. 作為一種分子力場模型,它們能夠以極高的準確性預測分子團簇中的二體和三體相互作用,這對于深入探討大分子系統的性質以及理解分子間相互作用的本質至關重要.

表3列出了MP2/aug-cc-pVTZ和本文方法計算的6個質子化精氨酸側鏈與水分子形成的團簇的分子間二體和三體相互作用能,以及本文方法與MP2/aug-cc-pVTZ的絕對誤差. 以二體作用能的計算結果為例,以MP2/aug-cc-pVTZ方法計算結果為基準,本文方法在計算這6個團簇的二體作用能時,絕對誤差的平均值為8.52 kJ·mol-1,相對誤差的平均值為 -1.22%. 這顯示出了本文方法在計算準確度方面的卓越表現. 具體而言,本文方法計算Cluster-B [(Arg(+) side chain· (H2O)10)]團簇的二體作用能與標準數據僅相差0.17 kJ·mol-1,而最大絕對誤差出現在Cluster-D[(Arg(+) side chain· (H2O)15)]團簇,其計算的二體作用能的絕對誤差為25.41 kJ·mol-1,相對誤差為 -3.9%.

表3 采用不同方法計算質子化精氨酸側鏈與水分子團簇的二體作用能、三體作用能及相應絕對誤差Table 3 Two-body interaction energies,three body interaction energies,and corresponding absolute errors calculated using different methods for protonated arginine side chain and water molecules clusters kJ·mol-1

不同于僅能描述分子間二體相互作用的非極化力場模型,本文方法基于誘導化學鍵偶極方案將極化作用引入模型. 采用簡單而有效的誘導偶極與其他偶極或電荷之間的相互作用勢來描述體系的三體及以上相互作用能.在三體作用能計算方面,本文方法也顯示出了較高的計算準確度,其計算結果的絕對誤差的平均值為 -9.36 kJ·mol-1. 這不僅顯示出了本方法能夠較好地描述分子間的三體作用,同時也和二體作用能計算誤差形成了有益的互補,能夠更準確地描述這類團簇分子間的總相互作用.

總的來說,本文方法在預測本文所研究團簇的二體和三體相互作用能方面表現出了高準確度,絕對誤差的絕對值最大為 9.87 kJ·mol-1,絕對誤差的平均值為 -0.83 kJ·mol-1. 這不僅反映了本文方法在分子間相互作用能計算方面的卓越表現,同時也表明了本文方法在三體相互作用能與二體相互作用能計算的誤差之間的互補性,能更準確地描述了這類團簇分子間的相互作用. 另外,值得一提的是,本文方法在不同研究系統中的計算誤差沒有明顯地隨體系規模增大而增加,這顯示出了它在應用和推廣方面的出色潛力.

為了更全面地評估本研究中新的誘導化學鍵偶極計算方法的可靠性,以及新的擬合參數在計算蛋白質分子中質子化精氨酸片段與水分子形成的團簇中的二體和三體相互作用能上的可靠性,重新構建了包含完整質子化精氨酸側鏈的水團簇分子結構.具體過程如下:首先,將質子化的精氨酸分子(兩端分別用氮甲基和乙?；舛孙柡?置于水盒子中,并使用NPT系綜下的GTH-PBE[9]方法進行了10 ps的分子動力學模擬;其次,在平衡后的模擬中,每1 ps保存一個結構,共獲得了5個質子化精氨酸分子與水分子形成的團簇. 需要強調的是,為了保證參考數據的計算準確度,僅截取距離質子化胍基C原子0.45 nm以內的水分子與質子化精氨酸殘基構成團簇. 這些團簇的結構如圖3所示.

圖3 質子化精氨酸二肽分子與水分子形成的團簇Fig.3 Clusters formed by protonated arginine dipeptide interacting with water molecules

分別采用RI-MP2/aug-cc-pVTZ、AMOEBAbio18[10-11]和本文方法計算上述5個團簇的總二體作用能和總三體作用能,計算結果列于表4中. 表中展示了作為基準的RI-MP2/aug-cc-pVTZ方法在計算這5個隨機團簇中二體和三體作用能的結果. 同時,本文方法和常用的AMOEBAbio18可極化力場方法的計算結果也列于表中. 此外,還列出了各組數據之間的絕對誤差,以便更清晰地評估這些計算方法的性能. 表中的Δ2b+3b為二體作用能和三體作用能絕對誤差之和.

表4 采用不同方法計算封端飽和的質子化精氨酸與水分子團簇的二體作用能、三體作用能及相應絕對誤差Table 4 Two-body interaction energies,three-body interaction energies and corresponding absolute errors calculated using different methods for protonated arginine dipeptide-water clusters kJ·mol-1

首先,對于Cluster-G和Cluster-H這兩個擁有相同水分子數量的團簇進行討論. 不論是使用AMOEBAbio18方法還是本文提出的方法,都能夠成功預測它們的相對穩定性. 不過,值得注意的是,相較于AMOEBAbio18方法,本文方法采用了一種更為高效的非迭代化學鍵偶極距方案,從而實現了相當高的計算準確度. 這也反映了本文方法在計算效率上的優勢. 此外,本文方法在對其他隨機團簇的二體和三體相互作用能進行計算時,同樣表現出了令人滿意的計算準確度,這也部分證實了該計算方法和參數的可靠性.

然而,需要特別注意的是,本文中使用了原有的基于中性氨基酸的主鏈參數來完成質子化精氨酸主鏈相關作用的計算. 這導致了在含有質子化精氨酸片段的水團簇的多體相互作用能計算中,相對于含有質子化精氨酸側鏈的水團簇,計算誤差較大. 為解決這一問題,我們計劃在未來的研究工作中進一步有針對性的優化精氨酸主鏈相關基團的參數,以期獲得更高的計算準確度. 這個優化過程將有助于更準確地描述蛋白質中精氨酸殘基與水分子的相互作用,從而提高模擬結果的可靠性.

總的來說,通過這些結果,不僅證實了本文方法在多體相互作用計算中的有效性,還展示了其在計算效率上的優勢. 然而,在應對其他氨基酸殘基的挑戰時,還將繼續努力,以進一步完善和改進本方法的計算準確度,使其更符合實際情況. 這將為分子模擬和計算化學領域提供更為可靠的工具和方法,有望在藥物設計和生物化學研究中產生廣泛的應用前景.

3 結 論

本研究通過引入一種新的誘導化學鍵偶極計算方法,成功應用于蛋白質分子中質子化精氨酸片段與水分子形成的團簇的相互作用能計算.結果表明,新方法在預測本文涉及的尺寸較大的含有質子化精氨酸側鏈、二肽的隨機水分子團簇的二體、三體總作用能的絕對誤差的平均值為 -8.3 kJ·mol-1,相對誤差的平均值為3.8%,表現出了良好的計算準確度.在計算量上還明顯小于常用的可極化力場方法.然而,在處理精氨酸主鏈相關作用時,也面臨一些挑戰,需要進一步改進相關參數以提高計算準確度. 這一研究為理解分子間相互作用,尤其是帶電荷片段間的相互作用提供了重要的工具和方法,對于藥物設計、生物物理化學等領域的研究具有積極的影響,為未來的科學研究提供有益的參考.