蘇敏儀，劉慧思，林海霞，王任小

1中國科學(xué)院上海有機(jī)化學(xué)研究所, 生命有機(jī)化學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 上海 200032

2中國科學(xué)院大學(xué)，北京 100049

3上海大學(xué)理學(xué)院化學(xué)系，上海 200444

1 引言

理解靶蛋白-藥物相互作用的基本規(guī)律是分子靶向藥物設(shè)計(jì)的基礎(chǔ)。靶蛋白-藥物相互作用可分別從熱力學(xué)和動(dòng)力學(xué)兩方面來表征。其中，結(jié)合熱力學(xué)常數(shù)(如平衡結(jié)合常數(shù)Ka)主要描述藥物分子與靶蛋白結(jié)合達(dá)到熱力學(xué)平衡時(shí)形成了多大比例的復(fù)合物。而結(jié)合動(dòng)力學(xué)常數(shù)(如結(jié)合速率常數(shù)kon、解離速率常數(shù)koff和駐留時(shí)間1,2)則描述了在給定時(shí)間內(nèi)復(fù)合物解離速率的快慢。親合性相同的復(fù)合物可能具有不同的結(jié)合速率和解離速率(見圖1)，因此動(dòng)力學(xué)參數(shù)可反映復(fù)合物之間不同的相互作用，可以為藥物分子設(shè)計(jì)提供更多有用的信息3。針對(duì)快結(jié)合和快解離的分子，其藥效持續(xù)時(shí)間較短，藥物占有率較低，則所需要的劑量也會(huì)越大。反之對(duì)于解離速度較慢的藥物分子，其在靶標(biāo)上的駐留時(shí)間較長；當(dāng)藥物濃度降低至飽和濃度以下時(shí)，其依然可以發(fā)揮藥效，可降低使用劑量；與此同時(shí)也會(huì)提高藥物的選擇性，減少副反應(yīng)的發(fā)生3,4。

目前大多數(shù)研究仍基于熱力學(xué)性質(zhì)(如平衡結(jié)合常數(shù)Ka)來衡量藥物與靶蛋白的相互作用強(qiáng)弱。用來獲得熱力學(xué)性質(zhì)的體外親合性實(shí)驗(yàn)(in vitrobinding assay)通常是在一個(gè)封閉系統(tǒng)中進(jìn)行，即藥物分子和靶蛋白的濃度在實(shí)驗(yàn)過程中不變1，這與生物體內(nèi)持續(xù)變化的生理環(huán)境并不相符。在一個(gè)開放的體內(nèi)系統(tǒng)(in vivo)中，一些慢結(jié)合或慢解離的藥物在常規(guī)實(shí)驗(yàn)時(shí)間內(nèi)難以達(dá)到平衡狀態(tài)。此時(shí)，測量藥物與靶蛋白的動(dòng)力學(xué)數(shù)據(jù)比熱力學(xué)數(shù)據(jù)更能體現(xiàn)藥物分子在體內(nèi)發(fā)揮的藥效4。近年來，越來越多的研究表明：藥物分子與靶標(biāo)結(jié)合的熱力學(xué)平衡性質(zhì)與藥物在體內(nèi)的藥理活性相關(guān)性較弱，而結(jié)合動(dòng)力學(xué)性質(zhì)與藥物在細(xì)胞和體內(nèi)的藥理活性則有很強(qiáng)的相關(guān)性3-6。例如，Heitman課題組研究了10個(gè)A2A腺苷受體激動(dòng)劑細(xì)胞內(nèi)功效與親合性和駐留時(shí)間的關(guān)系，發(fā)現(xiàn)該類激動(dòng)劑的細(xì)胞內(nèi)功效與駐留時(shí)間具有更好的相關(guān)性(R2=0.90)，其功效與抑制平衡常數(shù)Ki的相關(guān)性則很弱(R2= 0.13)5。因此，在現(xiàn)代藥物研發(fā)過程中，除了關(guān)注藥物與靶蛋白親合性的變化，研究者同時(shí)也需要關(guān)注藥物的結(jié)合動(dòng)力學(xué)參數(shù)，才能提高藥物研發(fā)的成功率。

圖1 藥物與靶蛋白結(jié)合過程中的自由能變化Fig.1 Free energy profile of the drug-target binding process.

由于測定藥物動(dòng)力學(xué)參數(shù)的實(shí)驗(yàn)方法有時(shí)間分辨要求，使得以動(dòng)力學(xué)性質(zhì)為導(dǎo)向的藥物分子設(shè)計(jì)需要發(fā)展新的實(shí)驗(yàn)方法，用來改善動(dòng)力學(xué)性質(zhì)測量結(jié)果的重復(fù)性和可靠性3。與此同時(shí)，研究者也迫切希望發(fā)展可以預(yù)測解離速率常數(shù)的計(jì)算方法，與實(shí)驗(yàn)方法相輔相成。目前常用的方法是采用消耗較多計(jì)算資源的動(dòng)力學(xué)模擬來研究藥物從靶蛋白上解離的途經(jīng)7。近年來，也有一些研究采用典型的定量結(jié)構(gòu)-活性關(guān)系(QSAR)策略來建立模型，并宣稱能獲得較好的預(yù)測結(jié)果8,9。本工作借助近年來方興未艾的機(jī)器學(xué)習(xí)方法，發(fā)展定量結(jié)構(gòu)-動(dòng)力學(xué)關(guān)系(QSKR)模型，實(shí)現(xiàn)預(yù)測藥物分子的解離速率常數(shù)，為基于結(jié)構(gòu)的藥物設(shè)計(jì)提供指導(dǎo)和幫助。

目前文獻(xiàn)中報(bào)道的此類QSKR模型并不多，主要圍繞某個(gè)特定靶標(biāo)開展研究工作。2016年重慶大學(xué)梅虎課題組發(fā)表了一項(xiàng)相關(guān)工作8。該工作通過采用3D分子力場Volsurf描述符及相關(guān)物理化學(xué)性質(zhì)對(duì)37個(gè)HIV-1的蛋白抑制劑結(jié)構(gòu)進(jìn)行特征描述，應(yīng)用偏最小二乘法(PLS)和支持向量機(jī)算法(SVM)構(gòu)建QSKR模型。該模型在測試集上預(yù)測解離速率常數(shù)的相關(guān)系數(shù)R為0.772。2018年Wade等人發(fā)表了一項(xiàng)類似工作，通過COMBINE分析方法及PLS構(gòu)建了一個(gè)QSKR模型9，其數(shù)據(jù)集包含了66個(gè)HSP90蛋白抑制劑。該模型在測試集上的相關(guān)系數(shù)R2為0.86。這些研究工作顯示：發(fā)展QSKR模型是一條可行的道路。但是在這些工作中用于構(gòu)建QSKR模型的數(shù)據(jù)集樣本量不夠，靶標(biāo)單一，限制了QSKR模型在藥物研發(fā)中的實(shí)際應(yīng)用范圍。

因此，發(fā)展更具通用性的QSKR模型是我們追求的目標(biāo)。為了實(shí)現(xiàn)這一目標(biāo)，在本工作中我們首先構(gòu)建了一個(gè)已知蛋白-配體解離動(dòng)力學(xué)實(shí)驗(yàn)數(shù)據(jù)的、規(guī)模足夠大且結(jié)構(gòu)多樣性好的數(shù)據(jù)集作為基礎(chǔ)，然后借助現(xiàn)代機(jī)器學(xué)習(xí)技術(shù)來構(gòu)建QSKR模型，實(shí)現(xiàn)基于蛋白-配體復(fù)合物結(jié)構(gòu)來預(yù)測配體分子的解離速率常數(shù)。最終得到的模型在結(jié)構(gòu)多樣性好的測試集上取得了良好的預(yù)測精度。

2 計(jì)算方法

本工作的主要流程如下：首先通過查閱文獻(xiàn)，收集有解離速率常數(shù)實(shí)驗(yàn)數(shù)據(jù)的配體及其對(duì)應(yīng)的蛋白-配體結(jié)構(gòu)信息。由于收錄到的大部分配體其復(fù)合物結(jié)構(gòu)未知，因此需要借助分子模擬技術(shù)預(yù)測它們與靶蛋白的合理結(jié)合模式。我們采用分子對(duì)接方法產(chǎn)生配體與靶蛋白的初始結(jié)合構(gòu)象，然后采用2 ns時(shí)長的非限制性動(dòng)力學(xué)模擬方法進(jìn)一步優(yōu)化蛋白-配體的相互作用模式。對(duì)于已知復(fù)合物晶體結(jié)構(gòu)的配體，我們對(duì)其復(fù)合物結(jié)構(gòu)采取2 ns時(shí)長的限制性動(dòng)力學(xué)模擬。獲得所有配體分子與對(duì)應(yīng)靶蛋白的結(jié)合構(gòu)象后，我們依據(jù)蛋白-配體復(fù)合物結(jié)構(gòu)描述符，采用隨機(jī)森林算法構(gòu)建了預(yù)測蛋白-配體動(dòng)力學(xué)解離速率的QSKR模型，在內(nèi)部以及外部測試集上取得了較好的預(yù)測效果。

2.1 數(shù)據(jù)收集

本課題組在過去十余年間，自主發(fā)展了蛋白-配體復(fù)合物數(shù)據(jù)庫PDBbind10，收錄了近2萬個(gè)已知親合性實(shí)驗(yàn)數(shù)據(jù)和晶體結(jié)構(gòu)的復(fù)合物。在維護(hù)PDBbind數(shù)據(jù)庫的過程中，我們積累了大量的相關(guān)文獻(xiàn)，可供收集蛋白-配體解離動(dòng)力學(xué)數(shù)據(jù)。我們首先采用關(guān)鍵字搜索方法(如kon、koff、off-rate、dissociation rate等)，從本組積累的3萬余篇文獻(xiàn)中篩選出2800余篇文獻(xiàn)。隨后通過人工閱讀文獻(xiàn)，收集文獻(xiàn)中報(bào)道的蛋白-配體解離速率常數(shù)koff及其實(shí)驗(yàn)手段和條件等相關(guān)信息，并記錄相應(yīng)的靶點(diǎn)信息(如PDB編號(hào)和UniProt ID等)。接下來，我們檢查并刪除數(shù)據(jù)集中的冗余樣本。如果同一復(fù)合物(即配體分子與靶蛋白分子均相同)存在多個(gè)koff數(shù)據(jù)，則優(yōu)先保留復(fù)合物晶體結(jié)構(gòu)已知的樣本，或優(yōu)先保留由Surface Plasmon Resonance技術(shù)測得數(shù)據(jù)的樣本。最終收錄了406個(gè)已知解離速率常數(shù)實(shí)驗(yàn)數(shù)據(jù)的配體分子以及對(duì)應(yīng)的靶蛋白信息。本工作采用CD-hit程序(v 4.5.4)11,12對(duì)數(shù)據(jù)集中的蛋白序列進(jìn)行聚類分析(相似性閾值設(shè)為90%)。聚類結(jié)果顯示該數(shù)據(jù)集樣本可分成114組，涵蓋了不同的蛋白類型，如腺苷受體A2A、MAPK14、嗜熱菌蛋白酶等，結(jié)構(gòu)多樣性豐富。圖2為該數(shù)據(jù)集中所有樣本的-lgkoff數(shù)值的分布圖。對(duì)該分布進(jìn)行Shapiro-Wilk正態(tài)分布檢驗(yàn)13,14，結(jié)果顯示該分布在95%置信水平下符合正態(tài)分布(p-value = 0.429)。

我們隨后采用ECFP4指紋15計(jì)算了該數(shù)據(jù)集中所有配體分子的結(jié)構(gòu)相似性矩陣。結(jié)果表明大多數(shù)配體分子之間的相似性非常低(見圖3)。同時(shí)，本工作計(jì)算了406個(gè)配體分子的一些基本物化性質(zhì)，包括分子量、可旋轉(zhuǎn)鍵數(shù)目、氫鍵受體數(shù)目、氫鍵供體數(shù)目、極性表面積和AlogP值。總體來看這些物化性質(zhì)的分布范圍較廣(見圖4)。以上所有結(jié)果提示該數(shù)據(jù)集具有顯著的多樣性。因此，采用該數(shù)據(jù)集為基礎(chǔ)最終導(dǎo)出的QSKR模型在原則上具有通用性。

圖2 原始數(shù)據(jù)集中406個(gè)配體分子的解離速率常數(shù)(koff)分布圖Fig.2 Distribution of the dissociation rate constants(koff) of the 406 ligands in the primary data set.

圖3 原始數(shù)據(jù)集中配體分子的相似性矩陣Fig.3 Similarity matrix of the ligand molecules in the primary data set.

圖4 原始數(shù)據(jù)集中406個(gè)配體分子的物化性質(zhì)分布Fig.4 Distribution of six physicochemical properties of the 406 ligands in the primary data set.

2.2 結(jié)構(gòu)處理

對(duì)于上述數(shù)據(jù)集中收錄的未知結(jié)合模式的配體分子，我們采用分子模擬方法預(yù)測其結(jié)合模式。首先對(duì)結(jié)構(gòu)進(jìn)行預(yù)處理：靶蛋白分子修補(bǔ)缺失殘基 采 用 Schr?dinger軟 件 (v.2018-3)的 Advanced Homology Modeling模塊，保存為PDB格式文件。預(yù)測配體分子質(zhì)子化狀態(tài)采用LigPrep模塊(假設(shè)pH= 7)，保存為Mol2格式。對(duì)于上述數(shù)據(jù)集中收錄的已知復(fù)合物晶體結(jié)構(gòu)的配體分子，采用ProToss工具對(duì)復(fù)合物結(jié)構(gòu)添加氫原子及預(yù)測質(zhì)子化狀態(tài)。ProToss是由Hamburg大學(xué)Rarey課題組2014年發(fā)表的軟件方法16。該軟件除了可以預(yù)測蛋白質(zhì)殘基的質(zhì)子化狀態(tài)及添加缺失的原子，還可根據(jù)配體周圍的殘基環(huán)境合理地預(yù)測配體分子的質(zhì)子化狀態(tài)、互變異構(gòu)體和氫原子坐標(biāo)。

2.3 預(yù)測配體與蛋白的結(jié)合模式

2.3.1 分子對(duì)接

針對(duì)那些未知結(jié)合模式的配體分子，我們采用分子對(duì)接技術(shù)獲得其與靶蛋白結(jié)合的初始結(jié)構(gòu)模型。我們假設(shè)具有相似化學(xué)結(jié)構(gòu)的配體分子如果與相同的靶蛋白結(jié)合會(huì)采用相似的結(jié)合模式。因此，可以采用結(jié)構(gòu)類似配體分子形成的復(fù)合物晶體結(jié)構(gòu)作為參考模板來開展分子對(duì)接。具體過程如下：首先采用Schr?dinger軟件中提供的Receptor Grid Generation工具，以相關(guān)晶體復(fù)合物中的配體為中心定義對(duì)接格點(diǎn)盒子。隨后選擇Ligand Docking工具，將未知結(jié)合模式的配體與其參考晶體蛋白進(jìn)行分子對(duì)接，采用GlidScore-SP17-19作為打分函數(shù)。完成分子對(duì)接后，我們對(duì)分子對(duì)接輸出的結(jié)果進(jìn)行檢查，挑選與參考模板結(jié)合模式相似且打分值較好的結(jié)合構(gòu)象，進(jìn)行后續(xù)的分子動(dòng)力學(xué)模擬步驟。如果分子對(duì)接輸出的高分值結(jié)合構(gòu)象與參考模板的結(jié)合模式不相似，則同時(shí)保留高分值以及與參考模板結(jié)合模式相似的對(duì)接構(gòu)象，進(jìn)行后續(xù)的分子動(dòng)力學(xué)模擬步驟。之后挑選由MM-GB/SA算法20計(jì)算的結(jié)合自由能更好的體系作為該配體分子與靶蛋白的合理結(jié)合模式。

2.3.2 分子動(dòng)力學(xué)模擬

在上述分子對(duì)接過程中，靶蛋白分子具有的柔性沒有考慮在內(nèi)。因此，我們采用模擬時(shí)間較短的分子動(dòng)力學(xué)模擬來進(jìn)一步優(yōu)化蛋白-配體復(fù)合物的分子對(duì)接模型。這里為了優(yōu)化晶體結(jié)構(gòu)可能存在的原子碰撞以及預(yù)測復(fù)合物體系中對(duì)配體解離過程起關(guān)鍵作用的水分子，蛋白-配體晶體結(jié)構(gòu)和由分子對(duì)接獲得的復(fù)合物結(jié)構(gòu)均進(jìn)行分子動(dòng)力學(xué)模擬：已知晶體結(jié)構(gòu)的復(fù)合物進(jìn)行2 ns時(shí)長的限制性動(dòng)力學(xué)模擬；分子對(duì)接預(yù)測的復(fù)合物結(jié)構(gòu)則進(jìn)行2 ns時(shí)長的非限制性動(dòng)力學(xué)模擬。

分子動(dòng)力學(xué)模擬采用AMBER軟件(2014版，GPU平臺(tái))21完成。首先使用Gaussian 09軟件22在HF/6-31G**水平計(jì)算配體分子的靜電勢(shì)分布，并采用限制性靜電勢(shì)擬合方法(RESP)23從靜電勢(shì)中擬合出配體的部分電荷。分子動(dòng)力學(xué)模擬中分別應(yīng)用FF14SB力場24和GAFF力場25來處理靶蛋白分子和配體分子。復(fù)合物置于正方體形狀的TIP3P水盒子26中，盒子表面到復(fù)合物結(jié)構(gòu)上任意原子的最小距離為10 ? (1 ? = 0.1 nm)，使用 Na+和Cl-中和整個(gè)體系的電荷。隨后聯(lián)合采用最陡下降法和共軛梯度算法分三步優(yōu)化體系：(1)約束蛋白-配體復(fù)合物，優(yōu)化水分子間的相互作用；(2)約束蛋白分子，優(yōu)化水分子和配體分子；(3)無約束條件下優(yōu)化整個(gè)體系。

接著將復(fù)合物體系從0 K加熱至300 K，加熱時(shí)間為200 ps；加熱結(jié)束后，在NPT條件模擬300 ps以平衡體系的密度。在加熱和平衡階段中，溶質(zhì)分子上的原子均施加一個(gè)數(shù)值為10 kcal·mol-1·?-2(1 cal = 4.1868 J)諧振約束，并采用SHAKE算法27對(duì)氫原子進(jìn)行限制。體系預(yù)熱平衡結(jié)束后執(zhí)行動(dòng)力學(xué)模擬。針對(duì)來源于復(fù)合物晶體結(jié)構(gòu)的體系，進(jìn)行2 ns時(shí)長的限制性動(dòng)力學(xué)模擬，蛋白分子上所有原子施加數(shù)值為10 kcal·mol-1·?-2的約束力；針對(duì)來源于分子對(duì)接的復(fù)合物體系，進(jìn)行2 ns時(shí)長的非限制性動(dòng)力學(xué)模擬(見圖5)。所有模擬中同樣采用SHAKE算法對(duì)氫原子進(jìn)行限制，模擬時(shí)間步長為2 fs。待分子動(dòng)力學(xué)模擬結(jié)束后，使用共軛梯度法對(duì)模擬最后時(shí)刻的結(jié)構(gòu)進(jìn)行優(yōu)化。最后使用AMBER軟件中的cpptraj模塊28提取優(yōu)化后的構(gòu)象，作為最終的復(fù)合物結(jié)構(gòu)用于后續(xù)步驟分析。有發(fā)表文獻(xiàn)指出水分子對(duì)配體的解離過程起著關(guān)鍵作用4。因此在提取蛋白-配體復(fù)合物構(gòu)象的同時(shí)，我們也保留其中的關(guān)鍵水分子，即氧原子同時(shí)與蛋白和配體的距離小于4 ?的水分子29。

由于本工作涉及多種類型不同的復(fù)合物體系，所以統(tǒng)一時(shí)長(2 ns)的動(dòng)力學(xué)模擬并不能確保所有體系達(dá)到相對(duì)穩(wěn)定的狀態(tài)。為此，我們采用cpptraj模塊進(jìn)行軌跡分析：如果某個(gè)復(fù)合物體系在2 ns時(shí)長內(nèi)未達(dá)到穩(wěn)定狀態(tài)，則該體系繼續(xù)進(jìn)行動(dòng)力學(xué)模擬直到相對(duì)穩(wěn)定為止。另外如前一小節(jié)所述，某些配體分子對(duì)應(yīng)多個(gè)可能的分子對(duì)接構(gòu)象，因此以上步驟會(huì)產(chǎn)生多個(gè)復(fù)合物結(jié)構(gòu)。在這種情況下，我們基于分子動(dòng)力學(xué)軌跡上最后1 ns進(jìn)行MM-GBSA計(jì)算，根據(jù)模擬結(jié)果和結(jié)合能計(jì)算結(jié)果挑選包含最合理結(jié)合模式的復(fù)合物結(jié)構(gòu)。

圖5 采用分子模擬得出兩個(gè)結(jié)構(gòu)相似的配體分子與同一靶蛋白的結(jié)合模式：(A)已知晶體結(jié)構(gòu)的復(fù)合物(PDB entry 3LDP)經(jīng)2 ns限制性動(dòng)力學(xué)模擬之后的結(jié)合模式；(B)未知復(fù)合物結(jié)構(gòu)的配體經(jīng)分子對(duì)接和2 ns非限制性動(dòng)力學(xué)模擬之后的結(jié)合模式Fig.5 The binding modes of two similar ligands derived by molecular simulation: (A) Binding mode after 2 ns restraint MD of a crystal structure (PDB entry 3LDP);(B) Binding mode after 2 ns unrestraint MD of the protein-ligand complex given by molecular docking.

2.4 應(yīng)用隨機(jī)森林算法構(gòu)建預(yù)測模型

通過以上步驟，我們?yōu)樵紨?shù)據(jù)集中所有406個(gè)蛋白-配體體系均構(gòu)建了對(duì)應(yīng)的復(fù)合物三維結(jié)構(gòu)模型。然后我們基于蛋白-配體相互作用結(jié)構(gòu)描述符，采用隨機(jī)森林算法來構(gòu)建預(yù)測配體分子解離速率常數(shù)的QSKR模型。

2.4.1 描述符的計(jì)算與選擇

在本工作中，我們采用蛋白-配體復(fù)合物結(jié)構(gòu)描述符來表征蛋白-配體相互作用。該套描述符源自Ballester等人的發(fā)表工作30。該套描述符統(tǒng)計(jì)配體分子與靶蛋白分子之間原子對(duì)的數(shù)目，考慮的蛋白原子類型包括C、N、O、S四種，配體原子類型包括C、N、O、P、S、F、Cl、Br、I九種。所以，理論上共存在36種配體-蛋白原子對(duì)，即基本描述符的總數(shù)為36個(gè)。對(duì)于每一種原子對(duì)(即每個(gè)基礎(chǔ)描述符)，則進(jìn)一步按照原子之間的距離來劃分。以蛋白分子上碳原子與配體分子上氧原子所形成的[C-O]原子對(duì)為例(圖6)：當(dāng)設(shè)定距離閾值為12 ?，分劃區(qū)間寬度為2 ?時(shí)，則分別統(tǒng)計(jì)原子間距離為[0, 2) ?、[2, 4) ?、[4, 6) ?、[6, 8) ?、[8, 10) ?、[10, 12) ?內(nèi)的[C-O]原子對(duì)的數(shù)目，共產(chǎn)生6個(gè)特征變量。由上所述，此類描述符可反映配體分子所處的環(huán)境，與在靶蛋白上的解離路徑相關(guān)。

圖6 示意圖：蛋白C原子(灰色)與配體O原子(紅色)形成原子對(duì)的統(tǒng)計(jì)方法Fig.6 Illustration of how the atom pairs formed between the carbon atoms (grey) on protein and an oxygen atom (red) on ligand are counted.

我們根據(jù)406個(gè)蛋白-配體復(fù)合物三維結(jié)構(gòu)模型產(chǎn)生描述符，并探索了不同距離閾值以及不同劃分區(qū)間寬度對(duì)最終QSKR模型預(yù)測性能的影響。具體情況見表1。當(dāng)距離閾值越大，劃分區(qū)間寬度越小，該套描述符可產(chǎn)生的特征變量則越多。而采用過多的特征變量會(huì)導(dǎo)致模型出現(xiàn)過擬合問題。因此，需要過濾其中識(shí)別能力較弱的特征變量，避免削弱模型的預(yù)測能力。我們采用方差法剔除無效或較弱的特征變量，過濾方差值小于設(shè)定閾值的變量。當(dāng)特征變量的方差很小時(shí)，即所有樣本在該特征變量的值變化較低，提示該特征變量的識(shí)別能力較差。本工作中我們分別考察了方差閾值分別為0、1和2時(shí)，所選擇的特征變量對(duì)最終QSKR模型預(yù)測能力的影響。

2.4.2 構(gòu)建QSKR模型

在本工作中，我們采用Python軟件提供的sklearn庫31(Anaconda發(fā)行平臺(tái)，版本4.3.0)中提供的隨機(jī)森林算法來構(gòu)建QSKR模型。該QSKR模型采用以上述描述符為自變量X，解離速率常數(shù)的負(fù)對(duì)數(shù)(-lgkoff)為因變量Y。評(píng)價(jià)QSKR模型性能的指標(biāo)包括：(1)衡量預(yù)測值與實(shí)驗(yàn)值之間相關(guān)性的Pearson相關(guān)系數(shù)R；(2)描述預(yù)測值與實(shí)驗(yàn)值相差程度的均方根偏差RMSE。此處需要說明的是：由于隨機(jī)森林算法過程中涉及使用隨機(jī)數(shù)，不同隨機(jī)數(shù)構(gòu)建的模型的預(yù)測性能表現(xiàn)具有差異性。因此，我們?cè)隍?yàn)證集和測試集上評(píng)價(jià)不同QSKR模型性能時(shí)，將構(gòu)建10個(gè)平行模型，取10個(gè)模型在測試集上評(píng)價(jià)結(jié)果的平均值來評(píng)判該QSKR模型的預(yù)測性能。

本工作中QSKR模型構(gòu)建流程如圖7所示。各主要步驟具體說明如下：

(1)數(shù)據(jù)集的分割：我們將數(shù)據(jù)集(樣本總數(shù)406)按3 : 1 : 1比例抽取樣本，分別組成訓(xùn)練集、內(nèi)部驗(yàn)證集和外部測試集。首先，隨機(jī)從整個(gè)數(shù)據(jù)集中抽取了81個(gè)樣本組成外部測試集。隨后，采用Kennard-Stone方法將剩余的樣本分割成訓(xùn)練集和驗(yàn)證集。Kennard-Stone算法是一種從數(shù)據(jù)集中挑選代表性子集的方法32,33，可實(shí)現(xiàn)均勻分割原始數(shù)據(jù)集描述符空間。該方法首先在數(shù)據(jù)集中找到歐幾里得距離最遠(yuǎn)的兩個(gè)樣本作為訓(xùn)練集成員。樣本間的歐幾里得距離即兩樣本間的描述符距離。接下來從剩余的樣本中挑選與前面所選樣本最近鄰距離最大者，并分配到訓(xùn)練集中。重復(fù)該步驟直至訓(xùn)練集的樣本量為244，剩余的81個(gè)樣本組成內(nèi)部驗(yàn)證集。

表1 計(jì)算蛋白-配體原子對(duì)描述符時(shí)所考慮的距離設(shè)置變量Table 1 Distance parameters considered in counting protein-ligand atom pair descriptors.

圖7 本工作中構(gòu)建QSKR模型的流程圖Fig.7 Workflow for deriving the QSKR model in our study.

(2)優(yōu)化超參數(shù)：我們考察了在不同距離閾值、劃分區(qū)間寬度及方差選擇水平下的描述符對(duì)模型預(yù)測性能的影響。在構(gòu)建模型時(shí)所定義的超參數(shù)會(huì)影響模型的性能表現(xiàn)。因此，基于每組描述符構(gòu)建的模型，我們首先采用5倍交叉驗(yàn)證優(yōu)化超參數(shù)，然后根據(jù)交叉驗(yàn)證結(jié)果挑選表現(xiàn)最好的超參數(shù)組合作為該組描述符下的超參數(shù)。

(3)模型的內(nèi)部驗(yàn)證：對(duì)于每組描述符采用步驟(2)中得到的超參數(shù)，在整個(gè)訓(xùn)練集中重新訓(xùn)練模型，并在內(nèi)部驗(yàn)證集上測試模型的預(yù)測性能。結(jié)合模型在驗(yàn)證集上的結(jié)果和在訓(xùn)練集中的交叉檢驗(yàn)結(jié)果，挑選表現(xiàn)最好的模型。

(4)模型的外部測試：采用測試集測試步驟(3)中得到的模型的預(yù)測性能。

這里需要說明的是：在本工作中我們事實(shí)上測試了sklearn庫中提供的4種常用機(jī)器學(xué)習(xí)算法，包括隨機(jī)森林、決策樹、支持向量機(jī)以及k-最近鄰算法。每種模型的構(gòu)建流程均采用上文中描述的基本步驟。結(jié)果顯示隨機(jī)森林算法構(gòu)建的模型在精度方面具有比較明顯的優(yōu)勢(shì)。因此本文中只報(bào)道隨機(jī)森林算法模型的相關(guān)結(jié)果，其他模型的結(jié)果則不再贅述。

3 計(jì)算結(jié)果與討論

由于我們建立的QSKR模型基于復(fù)合物的靜態(tài)結(jié)構(gòu)進(jìn)行計(jì)算，因此應(yīng)選擇一個(gè)在配體分子的解離過程中變化較小的描述符。與能量相關(guān)的描述符相比，本工作中所采用的原子對(duì)描述符具備這一特點(diǎn)29。在我們選擇的方法框架之下，為了獲得預(yù)測性能最優(yōu)的模型，我們分別探討了不同條件下產(chǎn)生、且經(jīng)不同特征選擇條件過濾后的描述符集對(duì)最終QSKR模型預(yù)測精度影響。

表2中總結(jié)了使用不同描述符集構(gòu)建的QSKR模型在驗(yàn)證集上獲得的預(yù)測值與實(shí)驗(yàn)值之間的相關(guān)性。從中我們可以看出：當(dāng)采用的特征來自距離閾值為15 ?、劃分區(qū)間寬度為3 ?、特征選擇方差水平為2時(shí)產(chǎn)生的描述符集構(gòu)建模型時(shí)，模型具有較好的預(yù)測性能，其在驗(yàn)證集上的相關(guān)系數(shù)為R= 0.671。這里我們將該模型記為模型1。獲得模型1的隨機(jī)森林算法中的超參數(shù)設(shè)置如下：隨機(jī)森林中回歸樹的數(shù)量為100；每棵回歸樹的特征數(shù)量為0.05；最大深度為60；使用out-of-bag策略采樣。在不同的特征方差水平下，來自距離閾值為15 ?、劃分區(qū)間寬度為3 ?的特征構(gòu)建的模型性能與其他采用不同距離條件構(gòu)建的模型相比，前者的性能表現(xiàn)較好。另一方面，來自距離閾值為15 ?，劃分區(qū)間寬度為3 ?的特征，在特征方差水平為2時(shí)所產(chǎn)生的模型較方差水平為0或1時(shí)所構(gòu)建的模型具有更優(yōu)異的預(yù)測性能，提示后者的描述集中存在識(shí)別能力較弱的描述符，削弱模型的預(yù)測性能。

隨后，我們對(duì)模型1進(jìn)行外部測試。模型1在外部測試集上的測試結(jié)果如表3所示。模型1表現(xiàn)出良好的預(yù)測能力，其預(yù)測值與實(shí)驗(yàn)值之間的相關(guān)系數(shù)R= 0.623，預(yù)測值均方根偏差RMSE = 1.06。我們還使用了另一個(gè)外部測試集來測試模型1。該數(shù)據(jù)集來自德國海森堡大學(xué)Wade課題組發(fā)表的工作，包含66個(gè)已知解離速率常數(shù)的HSP90蛋白抑制劑9。需要說明的是，該數(shù)據(jù)集有2個(gè)樣本與本工作的數(shù)據(jù)集重疊，故在我們的測試中對(duì)其進(jìn)行剔除，即該外部測試集由64個(gè)樣本組成。模型1在該外部測試集上仍然表現(xiàn)出了較強(qiáng)的預(yù)測能力(見表3)，其預(yù)測值與實(shí)驗(yàn)值之間的相關(guān)系數(shù)R=0.625，預(yù)測值均方根偏差RMSE = 0.98。該測試集只包含作用于單一靶點(diǎn)的配體分子，我們的QSKR模型在此測試集上表現(xiàn)出與多樣性測試集大致相當(dāng)?shù)木龋@是一個(gè)合理的結(jié)果。

表2 采用不同描述符集構(gòu)建的QSKR模型在驗(yàn)證集上的測試結(jié)果Table 2 Performance on the validation set by the QSKR model based on different descriptor sets.

表3 模型1和參照模型在不同數(shù)據(jù)集上的表現(xiàn)Table 3 Performance of model 1 and the reference model on different test sets.

為了進(jìn)一步評(píng)價(jià)我們得出的QSKR模型的真實(shí)能力，我們還構(gòu)建了一個(gè)簡單模型(Null Model)作為參照。該參照模型采用6個(gè)只與配體分子相關(guān)的物化性質(zhì)作為描述符，包括分子量、可旋轉(zhuǎn)鍵數(shù)目、氫鍵受體數(shù)目、氫鍵供體數(shù)目、極性表面積和AlogP值(見圖4)，采用相同的隨機(jī)森林算法(所有超參數(shù)均為軟件給定的缺省值)，在相同的訓(xùn)練集上訓(xùn)練得到。該模型在內(nèi)部驗(yàn)證集、外部測試集以及HSP90測試集上的表現(xiàn)請(qǐng)參見表3。結(jié)果顯示：我們構(gòu)建的QSKR模型的預(yù)測能力明顯強(qiáng)于該參照模型。這也說明單純用配體描述符不足以預(yù)測蛋白-配體復(fù)合物的解離速率常數(shù)，而本工作采用的蛋白-配體相互作用描述符則更適用于這一目標(biāo)。

4 結(jié)論

近年來，越來越多的研究表明藥物-靶標(biāo)結(jié)合的動(dòng)力學(xué)性質(zhì)與藥物在體內(nèi)的藥效有更強(qiáng)的相關(guān)性。發(fā)展能夠可靠預(yù)測藥物與靶蛋白結(jié)合動(dòng)力學(xué)性質(zhì)的計(jì)算模型，可以與實(shí)驗(yàn)方法相輔相成，在實(shí)踐中指導(dǎo)藥物分子的合理設(shè)計(jì)，從而提高藥物研發(fā)的成功率。

在本工作中，我們發(fā)展了一種通用型QSKR模型，可以基于蛋白-配體復(fù)合物結(jié)構(gòu)來預(yù)測配體分子的解離速率常數(shù)(koff)。我們首先從原始文獻(xiàn)中收集相關(guān)實(shí)驗(yàn)數(shù)據(jù)，構(gòu)建了一個(gè)含有406個(gè)樣本的數(shù)據(jù)集。該數(shù)據(jù)集的規(guī)模和多樣性都遠(yuǎn)遠(yuǎn)超過了前人工作中所采用的數(shù)據(jù)集。然后我們采用分子對(duì)接和分子動(dòng)力學(xué)模擬得出所有406個(gè)復(fù)合物的結(jié)構(gòu)模型，采用蛋白-配體原子對(duì)作為基本描述符，采用隨機(jī)森林算法構(gòu)建模型。在此過程中考察了多種距離閾值、劃分區(qū)間寬度以及特征選擇標(biāo)準(zhǔn)對(duì)模型性能的影響。綜合考慮在訓(xùn)練集以及驗(yàn)證集上得出的結(jié)果，最終確定當(dāng)距離閾值為15 ?，劃分區(qū)間寬度為3 ?，特征選擇方差水平為2時(shí)產(chǎn)生的描述符集合，其對(duì)應(yīng)的QSKR模型表現(xiàn)最好。在包含多樣性靶標(biāo)以及單一靶標(biāo)的外部測試集上，該模型均表現(xiàn)出良好的預(yù)測精度(R≈ 0.62)。此類QSKR模型雖不能在根本上解析配體分子解離的基本規(guī)律，但是在分子靶向藥物設(shè)計(jì)實(shí)踐中具有一定的應(yīng)用價(jià)值。另外，我們所得出的模型方法簡單，便于重復(fù)，有望為他人提供思路，促生性能更加卓越的計(jì)算預(yù)測方法。

Supporting Information:available free of chargeviathe internet at http://www.whxb.pku.edu.cn.