結合自由能微擾積分法及其在藥物研發中的應用*

甘雨滿，劉永宏,2，劉 鍇**

(1.廣西中醫藥大學海洋藥物研究院，廣西南寧 530200；2.中國科學院南海海洋研究所，中國科學院熱帶海洋生物資源與生態重點實驗室，廣東廣州 510301)

0 引言

分子動力學起源于1950年左右，其核心是運用牛頓(經典)力學模擬分子體系的運動，并從達到平衡的體系中抽取樣本，從統計學的角度計算體系的熱力學特征和其他宏觀性質[1,2]。2013年諾貝爾化學獎授予3位美國科學家：Martin Karplus、Michael Levitt 和 Arieh Warshel教授，以表彰他們發明的分子模擬方法對復雜生命科學領域發展的貢獻[3]。這是繼1998年諾貝爾化學獎頒給屬于理論計算的量子化學后，第二次將諾貝爾化學獎頒給理論計算，同時也標志著分子模擬方法得到實驗科學的認可。

分子動力學中最重要的是分子的力場(Force Field，FF)，它是計算分子各種性質的基礎，直接影響計算結果的可靠性及準確性。根據量子力學波恩-奧本海默近似，忽略電子的運動，可將分子體系的能量視為原子核位置的函數。與量子力學的從頭計算方法相比，分子動力學的計算量要小得多，因此能夠處理的計算體系也大得多；另一方面，在適當的范圍內，基于分子力場的計算精度與量子化學計算相差無幾。因此，對大分子復雜體系而言(如蛋白-小分子復合物)，基于分子力場的計算是一套經濟高效的計算方法。

眾所周知，藥物研發是一個耗資巨大且周期漫長的高風險過程。據估算，一個藥物成功上市需要約10年時間，花費高達14億美元[4,5]。近幾十年來，隨著分子生物學、結構生物學等學科的迅猛發展和高性能計算技術(軟件、硬件)的突飛猛進，分子動力學的MM/PB(GB)SA計算方法在先導化合物的篩選環節中得到廣泛應用[6]。MM/PB(GB)SA計算是利用分子動力學(Molecular Mechanism)并通過泊松-玻爾茲曼(Poisson Boltzmann)或廣義波恩(Generalized Born) 模型及分子表面積(Surface Area)的方法，近似求解靜電力(ELE)和范德華力(VDW)的貢獻。在蛋白-小分子復合物結合自由能的計算過程中，蛋白-小分子復合物的靜電力(ELE)和范德華力(VDW)的準確評估是計算的關鍵，其計算的具體步驟的相關綜述已非常全面詳細，在此不再累述[7,8]。

然而，鑒于MM/PB(GB)SA計算方法的局限性(動力學采樣、構型的熵變及介電常數選擇等)[6]，其計算精度及穩定性無法滿足藥物研發,尤其是藥物結構優化的高精度需求。例如：Sun等[9]系統研究1 508個蛋白-小分子復合物結構熵對MM/PB(GB)SA的影響，與實驗值相比，MM/PB(GB)SA計算的絕對結合自由能的平均偏差高達14.7 kcal/mol；Wang等[10]基于165個蛋白-小分子復合物的測試集計算，發現MM/PB(GB)SA計算與實驗值的平均絕對值偏差為9.7 kcal/mol。與之相反，結合自由能微擾積分法通過構建熱力學循環，可對蛋白-小分子復合物的結合自由能進行精確計算，越來越受到人們的重視。例如，He等[11]利用結合自由能微擾積分法計算134個蛋白-小分子復合物的結合自由能，其平均絕對誤差僅為0.71—0.94 kcal/mol，該精度可為先導化合物的結構優化提供非常重要的指導。本文首先介紹以分子動力學為基礎的結合自由能微擾積分法的計算過程，然后介紹限制勢的引入方法，最后通過具體的計算案例展現結合自由能微擾積分法計算在藥物研發中的應用前景。

1 結合自由能微擾積分法的計算過程

1.1 構建煉金術熱力學循環

通過分別模擬小分子在溶劑水和靶蛋白的口袋中，慢慢關閉或打開小分子與環境間的相互作用(ELE和VDW)，以上兩個過程構成了一個熱力學循環。由于小分子與環境間相互作用的關閉(或打開)的過程是不存在的非真實物理過程，因此該計算過程也被稱為“煉金術”。

圖1 煉金術熱力學循環示意圖[12]

1.2 結合自由能微擾積分過程

結合自由能微擾積分法，最早是由Kirkwood等[15]于1935年提出的。該方法通過向哈密頓量H(p,q)中引入一個耦合參數，再對進行積分計算出該過程的自由能變化，其中p和q分別代表離粒子的動量和原子位置。以圖1狀態A到B的過程中ELE變化為例，如果=1表示小分子處于狀態A，即小分子與環境具有完全的靜電力作用；λ=0代表狀態B，此時小分子與環境沒有靜電力作用，那么此過程中靜電力作用的能量變化為

在實際計算過程中，λ是從0到1間一系列連續的取值。每一個值，也被稱為一個窗口。首先創建一系列的λ點(或窗口)，分別計算每個窗口蛋白-小分子復合物的結合自由能；然后通過積分，求得此過程的自由能變化。如果λ的值越多(窗口越密集)，計算誤差將越小，但總體計算的耗時也將相應延長。一般來說，值的分布可采取等間距的辦法，如Aldeghi等[12]采用等間距dλ=0.1和0.05分別計算ELE和VDW的貢獻；λ值的分布也可以選擇非等間隔值，如He等[11]比較不同λ取值對計算的影響，并利用高斯求積法提出一個僅含9個λ的經濟計算策略。類似地，將靜電力換成范德華力時，也可求得狀態A到狀態B過程中范德華力的貢獻。最后，通過上述熱力學循環，小分子與蛋白的結合自由能可等價于從狀態A經狀態B—E等到達狀態F過程的自由能變化總和。

1.3 自由能計算精度的關鍵：動力學采樣

動力學模擬過程中的采樣是影響計算精度的一個關鍵因素[16,17]。以蛋白-小分子復合物為例，實驗上觀察到小分子進入蛋白結合口袋后，與蛋白穩定地結合在一起，其變化過程屬于自發行為。在能量上，其變化過程的吉布斯自由能為負值，根據范特霍夫等溫公式：

ΔG=-RTlnK，

其中R、T和K分別為氣體常數、反應的溫度和反應平衡常數。K代表了小分子與蛋白結合、解離兩種狀態間的動態平衡。

在實際的模擬過程中，由于逐漸關閉了小分子與蛋白及環境間的相互作用力(ELE和VDW)，此時蛋白與小分子復合物傾向于非結合狀態。如果要對蛋白-小分子的結合狀態有足夠的采樣，往往需要非常長的模擬時間。因此，為增強蛋白-小分子結合狀態時的模擬采樣，在實際計算過程中可通過人為引入限制勢來實現，即引入一個合理大小的作用力，約束蛋白-小分子復合物保持結合狀態。首先對蛋白與小分子復合物進行常溫常壓下(NPT)的常規動力學模擬，并對其平衡狀態進行一段時間的采樣(如10 ns)，評估蛋白-小分子復合物在平衡位置附近正常的波動范圍，然后分別從蛋白(P)和小分子(L)中分別選擇3個重原子(除氫原子以外的原子)作為參考原子建立限制勢。

如果原子允許的波動越小，則限制勢越強，對應地對復合物解離時的采樣時間越長，反之亦然。因此，在實際處理中僅考慮50%的原子的波動，試圖選取一個平衡點。首先，分別計算蛋白及小分子各重原子在平衡狀態中的正常波動(RMSF)，并按升序排列，選取前50%波動較小的原子建立參考列表。以小分子的3個參考原子選取為例，選擇RMSF最小的重原子為第一個原子LA，參考列表中距離LA最遠的原子為LB，參考列表中距離前兩個原子最遠的重原子為第三個參考原子LC。這樣做是確保對于柔性較大的分子，其質心也能保持穩定。與此類似地，從蛋白結構中選擇3個參考原子Pa、Pb和Pc(圖2)。限制勢則包含1個鍵長raA(即原子Pa和LA間的鍵長)、2個鍵角(θA、θB分別表示Pb-Pa-LA和Pa-LA-LB的鍵角)和3個二面角(φA、φB、φC分別表示二面角Pc-Pb-Pa-LA、Pb-Pa-LA-LB和Pa-LA-LB-LC)。

圖2 蛋白(P，灰色)-小分子(L，綠色)復合物限制勢中參考原子的示意圖

隨后，假定限制勢使得原子在平衡位置做簡諧振動。以兩個參考原子Pa和LA間鍵長(raA)為例，即簡諧勢U(x)=kr(x-x0)2，其中x表示任意一次采樣的鍵長raA，x0為兩個參考原子的平均平衡鍵長。那么簡諧常數kr為

其中，2σ表示上述采樣過程中，raA在室溫(T=298 K，其勢能為1kT，k為玻爾茲曼常數，即0.6 kcal/mol)下的正常波動范圍，從而可以求得鍵長的簡諧常數kr。依此類推，同樣可以求得其他參考原子間的鍵角、二面角間的簡諧常數(kθA、kθB、kφA、kφB、kφC)，從而可求得其限制勢(ΔF)大小[18,19]：

其中，k、T、V分別是玻爾茲曼常數、體系所處溫度和模擬體系的體積；raA.0、sinθA.0和sinθB.0分別是鍵長和鍵角處于平衡態時的平均值。

需要說明的是，增強動力學采樣是提高結合自由能計算精度的一個有效手段，但不是唯一方法。例如，Laury等[19]通過對葫蘆脲及14個不同大小化合物的結合自由能計算發現，使用可極化的分子力場也可以提高結合自由能微擾積分法的計算精度。

2 結合自由能微擾積分法在藥物研發中的應用

2.1 先導化合物結構的優化

先導化合物的結構優化是藥物能否進入臨床研究的關鍵。高精度的結合自由能微擾積分法可為先導化合物的結構優化提供重要的指導，加速藥物的研發進程。基于前述范特霍夫等溫公式可知，當蛋白-小分子復合物結合自由能相差0.6 kcal/mol時，蛋白-小分子復合物的反應平衡常數相差10倍。因此，高精度的結合自由能計算對精細地區分不同取代基的小分子與靶蛋白結合活性差異至關重要。表1選取了2015年以來結合自由能微擾積分法的一些代表性研究。當前主要關注結合自由能微擾積分法在計算案例中的計算精度，而總體上說，結合自由能微擾積分法的精度已經達到或接近于化學精度(1 kcal/mol)。如2015年，Wang等[20]針對8個靶蛋白系統地測試199個小分子的結合自由能，理論計算結果與實驗值的相關性達到0.75，且實驗誤差僅1.1 kcal/mol。該研究工作首次構建了1個高質量的訓練集，并含有多個靶蛋白及結構多樣的小分子化合物，為后續相關計算方法的研究提供了測試基礎與便利。另外，該研究也系統地評估分子對接(Glide 標準精度模式)、基于動力學的廣義玻恩模型(MM/GBSA)計算，以及結合自由能微擾積分法計算的穩定性和計算精度。與分子對接及MM/GBSA相比(分子對接和MM/GBSA計算結果與實驗的相關性分別為0.29和0.35)，結合自由能微擾積分法不僅精度高(0.75)，而且計算的魯棒性得以證實。Steinbrecher等[21]利用結合自由能微擾積分法計算實現對小分子碎片的活性預測，從而有利于早期先導藥物母核結構的開發。Aldeghi等[12]及Ciordia等[22]分別測試了蛋白-小分子復合物的結合自由能，其結果與實驗值誤差最小時僅差0.60和0.57 kcal/mol。Lenselink等[23]則將結合自由能的計算成功應用到G蛋白偶聯受體(GPCR)的藥物開發中，并獲得非常高的計算精度(均方根偏移RMSD僅為0.58—1.56)，為深入了解GPCR的結構和功能奠定了重要的理論基礎。Araki等[24]基于激酶蛋白構型的柔性，計算靶向周期蛋白依賴激酶2(CDK2)和細胞外調節蛋白激酶2(ERK2)的小分子抑制劑的結合自由能，并獲得非常好的相關性(R=0.81)。Li等[25]開發出一種基于高斯法則的結合自由能微擾積分法，對訓練集的計算表明其結合自由能變化與實驗活性的相關性R值達到0.69—0.94，并將該計算方法應用于靶向phosphodiesterase-10蛋白的苗頭化合物優化，最終得到優化后的化合物，并將其活性提高近2 000倍。

表1 結合自由能微擾積分法部分應用

續表1

2.2 蛋白-小分子復合物結合活性預測

結合自由能微擾積分法可對蛋白-小分子復合物結合活性進行理論預測，有利于篩選出最高活性的化合物，從而加速藥物研發效率。為此，Bhati等[26]通過考慮系統中的多樣本，提出一個提高積分法精度的計算方法，基于含有5個蛋白及55個小分子的測試集的計算結果顯示，其與實驗值的相關性R達到0.80—0.90。2017年Aldeghi等[27]通過結合自由能計算探討靶向Bromodomain (BRD)蛋白抑制劑的選擇性與活性，發現對化合物活性預測的平均誤差為0.81 kcal/mol，且相關性達到0.75；另外，該研究發現優化小分子磺酰胺基團的力場參數可提高計算準確性，也表明分子力場對于計算精度的重要性。Panel等[29]將結合自由能微擾積分法成功地應用到短肽-蛋白結合活性預測中，其計算的平均偏差最小達到0.37 kcal/mol，進一步拓展了結合自由能微擾積分法的應用范圍。

需要指出的是，盡管結合自由能微擾積分法展現出了高精度和穩定性，目前絕大部分的計算仍是基于回顧性分析的結果，在實際運用中的預測能力還需進一步觀察與評估。例如：Qian等[36]基于不同力場的結合自由能微擾積分法評估巨噬細胞游走抑制因子(Protein Macrophage Migration Inhibitory Factor)與其抑制劑間的結合自由能，計算過程中當小分子采用OPLS-AA/M力場及CM5原子電荷，基于蒙特卡洛或者分子動力學預測的結合自由能與實驗結果一致,分別為(-8.80±0.74)和(-8.46±0.85) kcal/mol，實驗值為(-8.98±0.28) kcal/mol；然而，當選用CHARMM和AMBER力場重復上述結合自由能計算時，預測的結合自由能最高誤差達6 kcal/mol。另一方面，Loeffler等[37]用結合自由能微擾積分法評估主流分子模擬軟件計算小分子與溶劑的溶解自由能(即結合自由能)的異同，發現盡管各軟件雖然在力場參數、數學表達式上不同，但對測試集的計算結果在數字上均符合預期，且其差別小于0.2 kcal/mol。以上結果說明，因不同的計算軟件和力場參數所導致的結合自由能預測計算誤差可能性非常小，6 kcal/mol的預測誤差很可能是由于蛋白與小分子的力場參數不匹配、不兼容造成的，另一方面也表明選擇匹配的力場參數對于準確計算至關重要。但如何才能選出適合的力場參數，卻仍是實際應用過程中的一個難點。盡管Aldeghi等[38]詳細介紹了結合自由能計算中的基礎方法與流程選擇，Abel等[39]提出改善結合自由能計算精度的一些經驗與策略，但目前還沒有一個統一的計算流程和方法可供不同軟件使用，仍需進一步發展與完善。因此，在實際運用中，需要對計算方法與流程、計算參數等進行測試與驗證，在此基礎上才能準確地預測蛋白-小分子復合物的結合活性。

2.3 脂水分配系數(logP)的計算

logP是藥物小分子一個非常重要的物理化學性質，對于藥物在體內的分布、代謝及吸收等至關重要。利用結合自由能微擾積分法對logP(或logD)的計算[40-44]，不僅對藥物的物理化學性質預測具有重要的應用價值，同時也有利于計算方法本身的完善。例如，logP計算已成為模擬計算測試(SAMPL)[44]中檢驗計算方法、力場參數的常規項目。當前，結合自由能微擾積分法對分子量較小的化合物的logP計算已經取得非常高的精度。例如：Garrido等[45]計算不同力場下的烷烴(n=1—8)的logP值，發現利用最優組合的參數條件(OPLS-AA/TraPPE力場)所得出的logP預測值與實驗值的平均誤差僅0.1個log單位。Bannan等[41]創建了一個含有41個小分子的測試集，并系統計算了該測試集在正辛烷及環己烷環境中的脂水分配系數，結果與實驗值的相關性(R)分別為0.70和0.83。

另外，對化合物logP的計算也可以為化合物小分子的結構改造與優化提供理論依據與合理解釋。例如，Liu等[46]利用結合自由能微擾積分法計算上市小分子藥物TAK-438[47,48](一個全新的鉀離子競爭型酸阻滯劑)不同氟取代位點的logP，發現logP與活性的變化存在良好的線性關系，認為單氟取代通過改變小分子化合物整體的偶極，從而調控化合物logP的變化，并成功地定性解釋了小分子抑制劑的不同氟取代點與抑制活性的關系。最后得出結論，單氟取代引起的logP變化是導致該系列化合物抑制活性巨大差異的根本原因。

2.4 蛋白與小分子結合模式的預測

蛋白與小分子結合模式是指小分子在靶蛋白的結合口袋中的活性構象及相互作用。一般來說，實驗室通過隨機生物學普篩或者盲篩，獲得少量具有生物活性的苗頭化合物；在此基礎上，利用X射線單晶衍射(X-ray)、核磁共振技術(NMR)或者冷凍電鏡(Cryo-EM)等技術手段解析出苗頭化合物與靶蛋白兩者間的相互作用模式。靶蛋白與苗頭化合物的結合模式是進一步對苗頭化合物進行合理地結構優化與修飾的基礎。但上述的實驗手段一般都具有實驗操作苛刻、耗時較長、花費較高等缺點,因此可以借助結合自由能微擾積分法快速準確地預測出蛋白與小分子的結合模式。鑒于結合自由能微擾積分法的密集型計算需求，可利用分子動力學對蛋白與小分子對接構型進行穩定性分析，快速剔除部分非穩定的對接構型。Liu等[49,50]從晶體結構數據庫中分別選出108和104個高質量的晶體結構作為測試集，并開展兩種方式的分子對接：第一種是自我對接，即將晶體結構中的小分子對接回原蛋白[49]；第二種是交叉對接，即用同一靶蛋白的不同抑制劑進行交叉對接[50]。針對每一個小分子分別選取3個代表性的對接構型(其中1個與晶體結構最接近的正確構型，2個誘餌構型)，通過動力學穩定性分析可剔除35%—55%的誘餌構型，從而大大降低下一步的蛋白-小分子結合活性預測的計算量。同時研究還發現，動力學穩定性分析還能夠優化分子對接的構型[49,50]，這一發現最近也被其他課題組所證實[51]。緊接著，Liu等[52]對第二個訓練集(共104個晶體結構)里剩余構象進行蛋白-小分子結合活性預測，發現72%(75/104)的復合物結合模式能夠被結合自由能微擾積分法正確地識別出。該結果表明利用分子動力學的結合自由能微擾積分法預測蛋白與小分子結合模式是可行的。最近，Zhang等[53]也成功將結合自由能微擾積分法應用于新型冠狀病毒(COVID-19)及其抑制劑瑞德西韋(Remdesivir)的結合模式預測的研究中，不僅預測了瑞德西韋可能的結合位點，同時也指出該結合模式下起重要作用的關鍵氨基酸，為進一步開發出更高活性的抑制劑提供理論指導和依據。另外，Aldeghi等[27]在預測靶向BRD蛋白家族的廣譜性小分子抑制劑的結合模式及其抑制活性時，發現預測活性與實驗值的誤差為1.76 kcal/mol，其相關性僅為0.48。該研究表明當前蛋白-小分子結合模式的預測難點，即分子動力學的結合自由能微擾積分法對蛋白-小分子兩者間的精細結構非常敏感，即便是相似的結合模式，不同蛋白-小分子結合活性預測值變化也比較大，因此需要對蛋白與小分子結合模式與活性預測仔細分析驗證。

3 展望

本文簡要介紹了以分子動力學為基礎的結合自由能微擾積分法在藥物研發中的應用。隨著高性能計算機的發展(如基于圖形處理器GPU的計算[11])，分子動力學模擬的計算速度也將迎來質的飛躍[54,55]。另一方面，越來越多高精度的量子化學的計算結果也為分子動力學力場參數的優化提供了高質量的數據來源[56,57]，而力場的不斷發展與完善也將為分子動力學的應用注入新的活力。

考慮到結合自由能微擾積分法計算步驟冗長、設置煩瑣，目前已有集成化、開源程序可供使用。例如：FEsetup[58]程序可兼容目前主流的動力學模擬軟件和程序，方便計算流程的自動化設置與計算，而alchemy-analysis[59]和MBAR[60]等分析工具則可簡化后期的數據處理與分析，如檢測動力學模擬過程中的采樣是否收斂、熱力學循環過程中的能量計算及誤差分析。這些程序與軟件的出現與發展，將大大降低結合自由能微擾積分法的使用門檻。相信隨著結合自由能微擾積分法計算的普及，其在藥物研發中將發揮更加重要的作用。