結(jié)合自由能微擾積分法及其在藥物研發(fā)中的應(yīng)用*

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

(1.廣西中醫(yī)藥大學(xué)海洋藥物研究院，廣西南寧 530200；2.中國(guó)科學(xué)院南海海洋研究所，中國(guó)科學(xué)院熱帶海洋生物資源與生態(tài)重點(diǎn)實(shí)驗(yàn)室，廣東廣州 510301)

0 引言

分子動(dòng)力學(xué)起源于1950年左右，其核心是運(yùn)用牛頓(經(jīng)典)力學(xué)模擬分子體系的運(yùn)動(dòng)，并從達(dá)到平衡的體系中抽取樣本，從統(tǒng)計(jì)學(xué)的角度計(jì)算體系的熱力學(xué)特征和其他宏觀性質(zhì)[1,2]。2013年諾貝爾化學(xué)獎(jiǎng)授予3位美國(guó)科學(xué)家：Martin Karplus、Michael Levitt 和 Arieh Warshel教授，以表彰他們發(fā)明的分子模擬方法對(duì)復(fù)雜生命科學(xué)領(lǐng)域發(fā)展的貢獻(xiàn)[3]。這是繼1998年諾貝爾化學(xué)獎(jiǎng)?lì)C給屬于理論計(jì)算的量子化學(xué)后，第二次將諾貝爾化學(xué)獎(jiǎng)?lì)C給理論計(jì)算，同時(shí)也標(biāo)志著分子模擬方法得到實(shí)驗(yàn)科學(xué)的認(rèn)可。

分子動(dòng)力學(xué)中最重要的是分子的力場(chǎng)(Force Field，F(xiàn)F)，它是計(jì)算分子各種性質(zhì)的基礎(chǔ)，直接影響計(jì)算結(jié)果的可靠性及準(zhǔn)確性。根據(jù)量子力學(xué)波恩-奧本海默近似，忽略電子的運(yùn)動(dòng)，可將分子體系的能量視為原子核位置的函數(shù)。與量子力學(xué)的從頭計(jì)算方法相比，分子動(dòng)力學(xué)的計(jì)算量要小得多，因此能夠處理的計(jì)算體系也大得多；另一方面，在適當(dāng)?shù)姆秶鷥?nèi)，基于分子力場(chǎng)的計(jì)算精度與量子化學(xué)計(jì)算相差無幾。因此，對(duì)大分子復(fù)雜體系而言(如蛋白-小分子復(fù)合物)，基于分子力場(chǎng)的計(jì)算是一套經(jīng)濟(jì)高效的計(jì)算方法。

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

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

1 結(jié)合自由能微擾積分法的計(jì)算過程

1.1 構(gòu)建煉金術(shù)熱力學(xué)循環(huán)

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

圖1 煉金術(shù)熱力學(xué)循環(huán)示意圖[12]

1.2 結(jié)合自由能微擾積分過程

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

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

1.3 自由能計(jì)算精度的關(guān)鍵：動(dòng)力學(xué)采樣

動(dòng)力學(xué)模擬過程中的采樣是影響計(jì)算精度的一個(gè)關(guān)鍵因素[16,17]。以蛋白-小分子復(fù)合物為例，實(shí)驗(yàn)上觀察到小分子進(jìn)入蛋白結(jié)合口袋后，與蛋白穩(wěn)定地結(jié)合在一起，其變化過程屬于自發(fā)行為。在能量上，其變化過程的吉布斯自由能為負(fù)值，根據(jù)范特霍夫等溫公式：

ΔG=-RTlnK，

其中R、T和K分別為氣體常數(shù)、反應(yīng)的溫度和反應(yīng)平衡常數(shù)。K代表了小分子與蛋白結(jié)合、解離兩種狀態(tài)間的動(dòng)態(tài)平衡。

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

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

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

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

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

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

需要說明的是，增強(qiáng)動(dòng)力學(xué)采樣是提高結(jié)合自由能計(jì)算精度的一個(gè)有效手段，但不是唯一方法。例如，Laury等[19]通過對(duì)葫蘆脲及14個(gè)不同大小化合物的結(jié)合自由能計(jì)算發(fā)現(xiàn)，使用可極化的分子力場(chǎng)也可以提高結(jié)合自由能微擾積分法的計(jì)算精度。

2 結(jié)合自由能微擾積分法在藥物研發(fā)中的應(yīng)用

2.1 先導(dǎo)化合物結(jié)構(gòu)的優(yōu)化

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

表1 結(jié)合自由能微擾積分法部分應(yīng)用

續(xù)表1

2.2 蛋白-小分子復(fù)合物結(jié)合活性預(yù)測(cè)

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

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

2.3 脂水分配系數(shù)(logP)的計(jì)算

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

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

2.4 蛋白與小分子結(jié)合模式的預(yù)測(cè)

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

3 展望

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

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