? ? ? ? ? ? ? ? ? ? ? ? ? ? gmx_MMPBSA的計算

? ? ? ? ?

? ? ? ? 前言：這一part將在閑時寫完。下一part將介紹線性相互作用能(LIE)方法，我感覺這個方法一直沒有受到重視，但是基于已有分子進行訓練來說，我覺得它可能也很靠譜。

? ? ? ? 背景

? ? ? ? ?gmx_MMPBSA是基于Amber中的MMPBSA.py進行的結(jié)合自由能計算，因此需要事先安裝AmberTool20或者更高版本。AmberTool20可由官方網(wǎng)站免費下載，大家可根據(jù)官方手冊說明自行安裝。它是基于AMBER MMPBSA.py用于計算終態(tài)自由能(gmx格式文件)的工具。包括兩個程序：gmx_MMPBSA：主程序，執(zhí)行計算功能。gmx_MMPBSA_ana：圖形界面分析數(shù)據(jù)并保存高質(zhì)量圖片。

? ? ? ? ? ? ? ??目前已知的支持力場的文件是：AMBER和CHARMM，所以在使用gromacs運算時，大家一定要選好力場。主推AMBER，GROMOS力場大家可以放棄了。

? ? ? ? ? 分子力場泊松玻爾茲曼表面積（Poisson?Boltzmann surface area，PBSA）和廣義波恩表面積（Generalized Born surface area，GBSA）是非常受歡迎的計算結(jié)合自由能的兩個模型。它們計算精確度高于對接函數(shù)，計算時間少于煉金術(shù)結(jié)合自由能（Alchemical binding free energy ABFE）計算。MM/PBSA和MM/GBSA被廣泛運用生物大分子的折疊、蛋白配體相互作用、蛋白-蛋白相互作用。?

? ? ? ? ? ? 熱力學上的自由能指的是體系內(nèi)部所具有的能量，決定了熱動力學過程的走向以及體系所處狀態(tài)的概率。研究自由能是研究生物分子運動和相互作用的手段之一，實驗測試自由能非常耗時且經(jīng)費投入較多，因此利用動力學模擬獲取分子體系的自由能變化是目前主流的計算方案，尤其是在藥物設(shè)計、蛋白質(zhì)折疊和蛋白-蛋白相互作用研究領(lǐng)域。在藥物與蛋白結(jié)合的熱力學過程中，在等溫等壓條件下體系的自由能通過如下公式定義：?F= -kBT lnZ。? 其中kB為玻爾茲曼常數(shù)，T為體系溫度，Z為體系分布函數(shù)，如果假定體系處于熱動力學穩(wěn)定狀態(tài)，則Z根據(jù)如下公式定義：

? ? ? ? ? ? 其中V0是體積常數(shù)，N是離子數(shù)，h是普朗克常數(shù)，P和V是壓力和體積。H(p,r)為指定動能（p）與位置（r）的粒子的漢密爾頓函數(shù)。有上述公式可知當原子數(shù)較少時，可以直接計算其絕對自由能，然而當原子數(shù)增高大時，計算絕對自由能變得非常困難。因此研究多原子體系的時候更多的是通過計算相對自由能變化來研究熱力學過程的走向。將研究狀態(tài)A與B之間的能量變化通過如下公式定義：

? ? ? ? ? ? ? 其中PA和PB為狀態(tài)A與B的概率。

? ? ? ? ? ? ? ?基于以上理論開發(fā)了很多計算藥物與蛋白結(jié)合自由能的計算機工具，對接函數(shù)速度快精度低，而更為復(fù)雜的煉金術(shù)結(jié)合自由能計算，例如自由能微擾（FEP）和熱力學積分算法（TI）精度高但計算非常耗時，而且由于生物分子體系的復(fù)雜度高，計算往往無法收斂，實際操作過程中另一個缺陷是準備步驟繁瑣，容易出錯。雖然近來有利用GPU加速來提升其計算速度，但是仍然不適合大規(guī)模篩選，只在具有較高計算資源項目組進行個別分子的優(yōu)化研究中得以運用。MM/PBSA和MM/GBSA方法則介于上述兩者之間，兼具精度速度，是廣泛被使用的計算結(jié)合自由能的工具。而且這些工具還能進行能量分解，每個氨基酸殘基對能量的貢獻值對于指導藥物設(shè)計也是極具幫助意義。

? ? ? ? ?MM/PB(GB)SA在藥物設(shè)計中的應(yīng)用

? ? ? ? ?1. 虛擬篩選

? ? ? ? 為了快速計算，通常采用簡單的打分函數(shù)來估計結(jié)合親和力，但是在大多數(shù)對接打分函數(shù)中簡化或者忽略了溶劑化自由能。在許多情況下，將分子對接和MM/PB(GB)SA結(jié)合是一種合理的策略，可正確識別結(jié)合構(gòu)象以及正確排序。有科研同仁同時對98個靶標的“對接能力”用MM/PB(GB)SA進行評估，發(fā)現(xiàn)在大多數(shù)情況下可以顯著提高找到配體的正確結(jié)合構(gòu)象的比例（成功率為69.4%）。但是預(yù)測的結(jié)果受到原子電荷、溶質(zhì)介電常數(shù)、MD模擬長度以及熵計算等等的影響。? ?

? ? ? ?2. 蛋白-蛋白相互作用

? ? ? ?蛋白-蛋白相互作用（PPIs）在活細胞的大多數(shù)生化過程中發(fā)揮關(guān)鍵作用。通過計算方法探索兩個蛋白間的相互作用，通常在初始階段使用精度低的打分函數(shù)快速篩選大分子數(shù)據(jù)庫，隨后使用更嚴格但更耗時的打分函數(shù)進行重新打分排序。研究表明，MM/GBSA比ZDOCK能更好地預(yù)測兩個蛋白結(jié)合的構(gòu)象。因此，考慮到低計算成本和相對較高的預(yù)測精度MM/GBSA可能是預(yù)測蛋白-蛋白系統(tǒng)正確結(jié)合構(gòu)象和結(jié)合親和力的強大工具。

? ? ? 3. 蛋白-核酸體系

? ? ? ??蛋白-核酸（DNA/RNA）相互作用在許多生物過程起到重要作用，如基因表達的調(diào)控、RNA剪切等。但是現(xiàn)有的蛋白-RNA相互作用（PRIs）打分函數(shù)并不可靠，因此準確預(yù)測PRIs的3D結(jié)構(gòu)和結(jié)合親和力仍然相當困難。MM/PB(GB)SA可被用于PRIs的研究。Orr等人使用MM/GBSA預(yù)測RNA修飾是否有利于與目標蛋白結(jié)合，預(yù)測結(jié)果與實驗數(shù)據(jù)有很高的相關(guān)性(r2> 0.9)。目前感覺不是很純熟，因此這方面的合作，一概不接。

問題集綿來啦，這都是重點：

[1] gmx_MMPBSA計算時幀數(shù)多少合適？

? ? ? ?在使用 gmx_MMPBSA 計算時，幀數(shù)的多少會影響計算結(jié)果的準確性和穩(wěn)定性。一般來說，幀數(shù)越多，結(jié)果越準確，但計算時間也會變長。因此，需要在時間和準確性之間取得一個平衡。通常建議使用至少1000個均勻分布的快照進行計算。這意味著從分子動力學模擬中選取至少1000個時間點，這些時間點應(yīng)該覆蓋整個模擬過程。如果模擬時間非常短，則可以使用更多的幀來覆蓋整個模擬過程。另外，如果您希望更加準確的結(jié)果，可以使用更多的幀來進行計算。如果計算時間不是主要問題，建議使用至少10,000個幀進行計算。如果需要更高的準確性，則可以使用更多的幀，但要記住，使用過多的幀可能會導致計算過于耗時。

? ? ? ? ?需要注意的是，使用過少的幀可能會導致結(jié)果不準確，而使用過多的幀可能會導致計算時間過長。因此，需要根據(jù)具體情況來確定合適的幀數(shù)。

? ? ? ? ?回答來自GPTchat，回答很中肯，采納。
? ? ? ? 后記：如何確定自己選了1000幀？

? ? ? ? ? ? ? ? ? ?答：一般我們步長是0.002ps，輸出精度我常選1000。這樣就是2ps保留一幀，1000ps=500幀，即1ns=500幀。我通常選后期10ns進行計算，即5000幀。在gmx_MMPBSA計算時選擇間隔5幀，即計算了1000幀。

[2] 關(guān)于多線程自己電腦實測：

? ? ? 普通MMGBSA和MMPBSA計算可以多線程。

? ? ? ?IE無法多線程。

? ? ? ?24核超線程，最多用12，選多報錯。

[3]GB和PB選擇哪個？

? ? ? ?Ⅰ：在應(yīng)用這些數(shù)據(jù)解析實驗現(xiàn)象或探究分子機制時，應(yīng)注意與實驗數(shù)據(jù)是否沖突，是否能解析得通，選擇更符合實際的。

? ? ? ?Ⅱ：?侯廷軍等人在2012年學術(shù)年會的一篇評估研究報告指出，MM/PBSA的預(yù)測結(jié)果更接近實驗數(shù)值，而MM/GBSA則對相對結(jié)合自由能預(yù)測性能更好。對MM/PB(GB)SA方法計算的結(jié)合自由能，應(yīng)結(jié)合實際進行取舍。侯廷軍, 李有勇. MM/PBSA和MM/GBSA對蛋白-配體自由能計算精度的評估研究. 中國化學會學術(shù)年會, 2012.??

? ? ? ?Ⅲ：MM/PB(GB)SA性能比較總結(jié)，選哪個自己決定：

? ? ? ? ?（1）MM/PB(GB)SA預(yù)測精度嚴重依賴于力場、電荷模型、連續(xù)溶劑模型、內(nèi)部介電常數(shù)、采樣方法以及構(gòu)象熵。? ??

? ? ? ? （2）力場在分子模擬中起到關(guān)鍵作用，決定了系統(tǒng)的所有相互作用。由于目前隱式溶劑模型都是非極化力場發(fā)展起來的，而MM/PBSA更適合極化力場。測試結(jié)果表明在MD模擬和MM/PB(GB)SA計算中使用相同的力場，不同力場的混合導致預(yù)測不準確。

? ? ? ? ?（3）結(jié)果表明PB和GB方法預(yù)測能力受到半徑的影響。

? ? ? ? （4）GB模型具有最好的預(yù)測結(jié)合親和力的能力。

? ? ? ? ? ? ? ? ? MM/PBSA方法比MM/GBSA方法表現(xiàn)出更好的排序性能。

? ? ? ? ? ? ? ? ? ?RESP電荷對MM/PB(GB)SA方法具有最佳性能，而基于AM1-BCC和ESP電荷預(yù)測的結(jié)果也相當令人滿意。

? ? ? ? ? （5）在計算絕對結(jié)合自由能的性能上，MM/GBSA預(yù)測能力比MM/PBSA差，但是MM/GBSA在不含金屬的體系中具有更好的結(jié)合親和力的排序能力。

? ? ? ? ? ?（6）個別體系中MMGBSA的精度已經(jīng)達到了FEP的計算精度，表明MMGBSA在計算結(jié)合自由能上的優(yōu)勢，相比之下MMPBSA的計算性能對體系敏感度更高，而且實際操作中MMPBSA的計算更加耗時，由于溶劑模型的影響MMPBSA經(jīng)常出現(xiàn)結(jié)合自由能為正數(shù)的異常情況，也是導致MMPBSA使用不如MMGBSA頻繁的原因之一（網(wǎng)上常問問題，原因就是這個，你算正無外乎選擇錯分組或者換MMGBSA）。

? ? ? （7）MD時長不是越長越好，對于MMGBSA計算1ns的動力學模擬很多時候也可以滿足計算需求。所以說還是有時候。

? ? ?（8）PB的原理是最嚴格的，所以大量GB模型的文章都用PB的結(jié)果當參照判斷GB的設(shè)置、求解有效Born半徑的方法的合理性。-來自盧天。

? ? ? ? ?參考文獻：[1] Wang E, Sun H, Wang J, et al. End-point binding free energy calculation with MM/PBSA and MM/GBSA: strategies and applications in drug design. Chemical Reviews, 2019, 119(16): 9478-9508.??