談談分子動力學模擬蛋白-配體復合物過程中配體發生脫離的原因

談談分子動力學模擬蛋白-配體復合物過程中配體發生脫離的原因

文/Sobereva@北京科音   2022-Jan-13

有不少人問怎么他做基于經典力場的蛋白質-配體復合物的分子動力學模擬的過程中配體跑出去了。鑒于這個問題被問得很頻繁，筆者遂寫個小文把常見的可能因素說明一下，便于讀者排查原因。

首先要憑基本化學常識和直覺判斷配體是否真有可能穩定與蛋白結合，可以在VMD程序中把配體周圍一定距離的蛋白質原子都顯示出來觀察，這需要利用within語句來定義顯示范圍，詳見《VMD里原子選擇語句的語法和例子》（http://www.shanxitv.org/504）。配體和蛋白結合常見作用機制有氫鍵（通常主要本質是靜電吸引作用，見http://www.shanxitv.org/513）、鹽橋（配體顯離子性部分與氨基酸側鏈帶相反電荷基團的靜電吸引）、范德華作用（具體來說是指色散吸引效應，pi-pi堆積本質上與此也相同）、疏水作用（本質上來自于溶劑的熵效應，不懂的話Google）。如果你對弱相互作用類型不太了解，可以看看此文里相關科普：《談談“計算時是否需要加DFT-D3色散校正？”》（http://www.shanxitv.org/413）。基本上，對弱相互作用有一定了解、接觸過較多蛋白-配體復合物體系研究的人都能憑感覺初步估計在一個給定的復合物結構中，一個小分子和蛋白質是否有可能穩定結合上。如果在初始結構下，上述這些導致配體和蛋白相互結合的作用非常不顯著，顯然小分子可能由于熱運動跑到溶劑里去，或者跑到蛋白質上其它更容易結合當前配體的區域去。值得一提的是，有一些程序可以幫助你大致了解一下配體和蛋白質之間的相互作用，比如LigPlot+程序（https://www.ebi.ac.uk/thornton-srv/software/LigPlus/）和https://proteins.plus在線服務器里的PoseView。此外，proteins.plus里的DoGSiteScorer還可以找出蛋白質中的口袋并圖形化展現出來，并給出Drug Score體現口袋的類藥性（結合藥物分子的可能性），有時能幫助你做一些粗略的判斷、輔助確定對接需要考慮的空間范圍。

初始的復合物結構是怎么得到的，顯然直接關系到跑動力學期間配體能否在初始位置附近呆住。有這么幾個情況，應當分情況考慮
(1)用的是X光衍射測的復合物結構：這樣的初始結構通常質量較高，小分子理應能在動力學過程中結合在初始的位置。但是也要注意，晶體狀態的復合物結構只是實驗測定條件下的最穩定結構，而在現實條件的動力學模擬過程中，由于熱運動和溶劑效應，小分子未必就一定能始終穩定呆在晶體中原本的位置。另外，有些晶體結構的解析度較低，比如三點幾埃，此時配體、氨基酸殘基位置的不確定性很高，可能此結構下側鏈與配體的相互作用情況和實際有明顯偏差，進而可能導致模擬開始后沒多久配體就跑掉。
(2)用的是分子對接得到的復合物結構：如果沒有復合物晶體結構但是有蛋白質晶體結構，通常就是用對接得到的打分最高到的復合物結構當初始結構。受制于采樣充分程度、打分函數的準確度、對受體的柔性考慮等因素，對接產生的初始結構的可靠性明顯不如較高解析度的晶體衍射測出來的結構。很多人做過分子對接后都是再專門跑一下分子動力學檢驗配體能否基本維持在對接產生的位置，從而檢驗對接產生的結構的合理性。使用這種結構當初始結構時，如果配體在動力學模擬過程中跑掉了，很可能就是因為對接產生的打分最高的結構并不合理。此時可以考慮用你覺得看上去也靠譜但打分不是最高的復合物結構試試，或者換其它對接程序得到的結構試試。
(3)用的是對接得到的結構，而且蛋白質結構是預測出來的：這無疑是可靠性最低的結構，因為蛋白質本身結構的可靠性都成問題。這種情況配體如果在模擬中跑出來，除了考慮對接的問題外，當然也得考慮蛋白質三維模型是否真的合理。蛋白質預測程序/方法很多，其中一種不理想時可以再考慮別的，選擇余地較大，比如基于同源模建的Modeller、Robetta等，基于深度學習預測蛋白質結構的AlphaFold 2和RoseTTAFold。
用上面(2)、(3)方式得到的初始結構，在跑動力學之前記得仔細觀察配體和受體結合區域的特征，憑常識估計是否真有可能靠非共價的方式穩定結合，明顯沒什么希望的結構直接pass掉。

用的力場不合理也是可能導致小分子在動力學過程中跑掉的原因，但記住這絕對不是唯一原因，別一發現小分子跑掉就總是賴力場。目前模擬蛋白質+小分子體系的常見組合是AMBER（蛋白質）+GAFF（小分子）、CHARMM（蛋白質）+CGenFF（小分子）、GROMOS（直接描述蛋白質，并結合ATB創建的配體拓撲文件）、OPLS-AA（直接描述蛋白質，并結合LigParGen創建的配體拓撲文件）。其實以上這些蛋白質+小分子描述的組合，對于前述的范德華作用的描述沒太大差異（而且通常范德華作用對蛋白質-配體結合的貢獻相對來說較次要），而且也不直接影響疏水作用（疏水作用是在動力學模擬過程中自然而然體現的，不是在力場層面上描述的），也沒法說哪個一定更好。對于GROMACS用戶，我個人建議優先考慮AMBER14SB結合acpype創建的基于GAFF力場的小分子拓撲文件，主要是因為這倆力場都比較魯棒，AMBER14SB有GROMACS現成的力場包而且描述蛋白質構象較好（雖然目前也有更新的AMBER19SB，但截止到撰文時網上還沒有現成的GROMACS力場包，對于描述蛋白-配體作用方面它比AMBER14SB沒有額外優勢），而且acpype程序用起來方便，還有在線版。對小分子拓撲文件創建不了解的話看《幾種生成有機分子GROMACS拓撲文件的工具》（http://www.shanxitv.org/266）。靜電作用對于蛋白-配體結合的穩定程度起至關重要的影響，而且它直接決定了對氫鍵和鹽橋的描述合理性。以上力場中描述靜電作用都是靠原子電荷基于庫侖公式來計算的，而這些力場對蛋白質的標準氨基酸殘基都已經明確定義了原子電荷，因此在特定的蛋白&配體力場組合下，影響結合穩定性的關鍵在于配體的原子電荷的確定，這是非常有講究的。對于AMBER+GAFF的組合，給配體用RESP2(0.5)原子電荷即便不是最理想至少也是最理想之一，可以用Multiwfn程序（http://www.shanxitv.org/multiwfn）基于Gaussian或ORCA等量子化學程序優化的結構和產生的波函數直接算出來，見《RESP2原子電荷的思想以及在Multiwfn中的計算》（http://www.shanxitv.org/531）。RESP2(0.5)對于其它蛋白-配體力場組合不一定是最理想的，因為牽扯到與范德華參數、二面角參數的耦合等問題，這里不打算展開細說，但至少用RESP2(0.5)仍是合理選擇之一。使用上述力場組合（假設蛋白質力場部分都是較新版本），對于蛋白質結構的描述都沒有什么問題，但對于成鍵方式不夠典型的配體來說，有可能出現在某個力場下它的結構能正確維持，而在另一個力場下則明顯不合理的情況。如果配體自身結構都沒法被正確描述（如某個局部本應是平面的而在動力學過程中卻嚴重彎曲了），自然也就不能指望配體和蛋白口袋的相互作用肯定能合理描述。因此，做動力學的時候要觀察軌跡，留意小分子在模擬過程中的結構是否和期望一致，必須確保沒有明顯異常。如果你沒有化學直覺，不知道一個比較特殊的配體分子（假設是有機的）理應是什么樣，可以用比如B3LYP-D3(BJ)/6-311G*這種常見也靠譜的幾何優化級別用量子化學程序做優化，還可以用molclus（http://www.keinsci.com/research/molclus.html）做構象搜索。但也應注意，對于柔性較大的配體，蛋白質與它的相互作用往往會明顯影響其構象，可能在結合時它的構象明顯和孤立狀態下自由能最低構象不同，這是很正常的。

如果在AMBER14SB+GAFF結合RESP2(0.5)原子電荷下模擬蛋白-配體復合物時出現配體中途跑掉的情況（而且看上去小分子自身結構的描述沒啥問題），換力場也大概率解決不了。雖然也可以再試一個別的力場組合，但如果發現還不行，一般也沒必要再試了，而應該從本文提及的其它角度考慮跑掉的原因。

有些特殊的蛋白-配體相互作用不是靠以上常見的力場能直接描述的。比如碘代苯和蛋白質上負電性原子間可以形成顯著的鹵鍵，這光靠原子電荷是沒法展現的，而必須表現出鹵原子周圍電荷分布的各向異性，這需要修改原先的模型，考慮額外點電荷，或考慮原子多極矩。例如在JCTC, 14, 5383 (2018)有人提出給GROMOS力場的鹵原子的σ-hole區域增加額外的正電性的點電荷使得鹵鍵能夠定性正確表現出來。PS：在GROMACS中如果想用此文的方式，需要依靠虛擬點來實現。另外，如果蛋白-配體之間有涉及過渡金屬的配位鍵出現，光靠原子電荷描述的靜電吸引作用一般也不足矣正確描述和維持之，此時往往需要加上限制勢人為強行維持住。而對于涉及到高價主族金屬離子的蛋白-配體結合，光靠原子電荷表現靜電作用可能也描述不充分，因為明顯的極化作用對結合會產生穩定化，這一般需要量子化學方法才能準確描述，因此這種情況下若配體在動力學過程中跑飛了的話也得考慮加上限制勢來強行維持住金屬離子與周圍原子的作用。

配體、殘基側鏈在結合時所處的質子化狀態也是要恰當考慮的，否則無法正確表現蛋白-配體間的實際相互作用。有的配體、氨基酸殘基在不同pH下有不同的質子化狀態，而且蛋白-配體的相互作用還可能影響質子化狀態。對于組氨酸尤其要注意考慮，它在模型體系中的pKa接近7，在生理環境下模擬時其側鏈的質子化態有較大不確定性，受環境影響明顯，既可以處于中性（側鏈上的一個氮被質子化），也可以帶一個正電荷（側鏈上兩個氮都被質子化）。當配體與組氨酸可能存在氫鍵作用時，組氨酸的質子化狀態應當是有利于氫鍵形成的。你缺乏相關知識和判斷能力的話也不要緊，有一些程序能自動給你判斷蛋白和配體合理的質子化狀態，例如https://proteins.plus在線服務器里的Protoss工具（介紹見https://www.zbh.uni-hamburg.de/forschung/amd/software/protoss.html）。如果你是用GROMACS的話，pdb2gmx產生的蛋白拓撲文件時可以加上-his選項來人工選擇各個組氨酸的質子化態成為Protoss給出的情況，而使用小分子拓撲文件產生工具時應當提供合理的質子化態的結構文件來產生。

動力學模擬設置不合理也是可能導致小分子跑掉的原因。做蛋白-配體的正式動力學模擬之前一般都是在給蛋白質和配體加上位置限制勢的情況下跑一段短暫的動力學，以使得溶劑分子在這個過程中能充分弛豫、迎合復合物結構。如果沒做這一步而直接跑常規的動力學的話，有可能會由于初始加的水的位置不夠理想，不僅沒有描述充分實際的水環境，甚至造成出現一些不合理的互斥，導致配體跑出去。另外，做生物體系的動力學我一般建議讓參考溫度從0 K開始緩慢線性上升，這樣最為穩妥。如果讓程序隨機產生初速度，而且初速度對應的溫度又不太低且沒加限制勢，也沒準會導致在模擬初期配體和蛋白的相互作用尚未充分穩定的時候配體由于擁有過高的動能而跑走。

如果復合物結構是通過對接方式產生的，發現動力學模擬還沒有跑多久配體就脫離了，可以考慮給小分子和與之相互作用的殘基之間恰當地加上一個或多個較弱的距離限制勢（怎么加根據具體結構隨機應變，比如如果配體和小分子之間通過氫鍵維持，可以給氫鍵作用的原子間加上限制勢維持住氫鍵距離），避免小分子跑掉。等動力學跑一陣子后，蛋白質的口袋區域的構象就弛豫了，蛋白-配體相互作用會比一開始的狀態有所增強，之后再去掉限制勢進行正式的模擬，說不定配體就能一直呆住了。注意，除非是由于前述的力場的局限性而不得不加限制勢，否則正式模擬時限制勢必須去掉，否則相當于給模擬引入了明顯的人為因素，發文章的話容易被審稿人抨擊。

不要以為配體在現實中就應該始終和蛋白是維持結合狀態的。配體L和蛋白P結合成復合物LP的速率為k_on[L][P]，其中k_on為結合速率常數；而LP解離為L和P的速率為k_off[LP]，其中k_off為解離速率常數。當達到穩態時滿足k_on[L][P]=k_off[LP]，并且有Kb=k_on/k_off=[LP]/([L][P])=1/Kd的關系，其中Kb稱為結合常數，Kd稱為解離常數。顯然，復合物的形成和復合物的解離是同時進行的，不是說配體結合到蛋白后就一直不出來了，也因此Kb不是個無窮大的值。k_off越大體現配體越容易解離、能結合的平均時間越短，這和解離過程的自由能壘有密切聯系。k_off較大的配體說明它結合得就是不很穩定，就是容易在動力學過程中跑出去，所以不能一味地認為模擬中出現解離的現象就一定是錯誤的。有的人跑復合物的模擬并非才跑了沒多久配體就脫離了，而是跑了好幾百ns才出現配體脫離現象，很可能正是體現出真實的配體解離現象。

還應當考慮自己做的模擬是否有真的意義，期望出現的現象是否真的是真實的。比如之前我看到有人研究的是酶，問怎么小分子在里面沒法一直呆住。實際上要求小分子呆住本來也沒意義，本身底物小分子偶然進去之后被催化完就走掉了，何故期望小分子能穩定呆在里面？