應用機器學習方法構建藥物分子解離速率常數的預測模型

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

1中國科學院上海有機化學研究所, 生命有機化學國家重點實驗室, 上海 200032

2中國科學院大學，北京 100049

3上海大學理學院化學系，上海 200444

1 引言

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

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

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

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

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

因此，發展更具通用性的QSKR模型是我們追求的目標。為了實現這一目標，在本工作中我們首先構建了一個已知蛋白-配體解離動力學實驗數據的、規模足夠大且結構多樣性好的數據集作為基礎，然后借助現代機器學習技術來構建QSKR模型，實現基于蛋白-配體復合物結構來預測配體分子的解離速率常數。最終得到的模型在結構多樣性好的測試集上取得了良好的預測精度。

2 計算方法

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

2.1 數據收集

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

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

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

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

圖4 原始數據集中406個配體分子的物化性質分布Fig.4 Distribution of six physicochemical properties of the 406 ligands in the primary data set.

2.2 結構處理

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

2.3 預測配體與蛋白的結合模式

2.3.1 分子對接

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

2.3.2 分子動力學模擬

在上述分子對接過程中，靶蛋白分子具有的柔性沒有考慮在內。因此，我們采用模擬時間較短的分子動力學模擬來進一步優化蛋白-配體復合物的分子對接模型。這里為了優化晶體結構可能存在的原子碰撞以及預測復合物體系中對配體解離過程起關鍵作用的水分子，蛋白-配體晶體結構和由分子對接獲得的復合物結構均進行分子動力學模擬：已知晶體結構的復合物進行2 ns時長的限制性動力學模擬；分子對接預測的復合物結構則進行2 ns時長的非限制性動力學模擬。

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

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

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

圖5 采用分子模擬得出兩個結構相似的配體分子與同一靶蛋白的結合模式：(A)已知晶體結構的復合物(PDB entry 3LDP)經2 ns限制性動力學模擬之后的結合模式；(B)未知復合物結構的配體經分子對接和2 ns非限制性動力學模擬之后的結合模式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 應用隨機森林算法構建預測模型

通過以上步驟，我們為原始數據集中所有406個蛋白-配體體系均構建了對應的復合物三維結構模型。然后我們基于蛋白-配體相互作用結構描述符，采用隨機森林算法來構建預測配體分子解離速率常數的QSKR模型。

2.4.1 描述符的計算與選擇

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

圖6 示意圖：蛋白C原子(灰色)與配體O原子(紅色)形成原子對的統計方法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.

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

2.4.2 構建QSKR模型

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

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

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

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

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

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

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

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

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

3 計算結果與討論

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

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

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

表2 采用不同描述符集構建的QSKR模型在驗證集上的測試結果Table 2 Performance on the validation set by the QSKR model based on different descriptor sets.

表3 模型1和參照模型在不同數據集上的表現Table 3 Performance of model 1 and the reference model on different test sets.

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

4 結論

近年來，越來越多的研究表明藥物-靶標結合的動力學性質與藥物在體內的藥效有更強的相關性。發展能夠可靠預測藥物與靶蛋白結合動力學性質的計算模型，可以與實驗方法相輔相成，在實踐中指導藥物分子的合理設計，從而提高藥物研發的成功率。

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

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