林開東 林曉倩2) 林緒波?

1) (北京航空航天大學(xué),醫(yī)學(xué)科學(xué)與工程學(xué)院/生物與醫(yī)學(xué)工程學(xué)院,北京市生物醫(yī)學(xué)工程高精尖創(chuàng)新中心,北京 100191)

2) (北京航空航天大學(xué)沈元學(xué)院,北京 100191)

1 引言

阻斷免疫檢查點蛋白與其配體的結(jié)合作為腫瘤治療的方法之一,近年來在臨床應(yīng)用中迅速發(fā)展.正常生理狀態(tài)下,部分負(fù)調(diào)節(jié)因子作為免疫檢查點來抑制T細(xì)胞的過度激活,確保免疫反應(yīng)保持自我耐受[1].然而不幸的是,腫瘤細(xì)胞可以利用這種機制誘導(dǎo)T細(xì)胞衰竭,形成促進(jìn)腫瘤生長和侵襲的微環(huán)境,從而逃避免疫系統(tǒng)的攻擊[2-4].為了重新激活并增強T細(xì)胞介導(dǎo)的抗腫瘤功能,前人已經(jīng)設(shè)計了一系列療法來阻斷免疫檢查點蛋白與其配體的結(jié)合[5,6].其中,細(xì)胞程序性死亡蛋白I(programmed cell death protein 1,PD-1)是最受關(guān)注的免疫檢查點之一,其通常表達(dá)于活化的T細(xì)胞、自然殺傷性細(xì)胞、B淋巴細(xì)胞和其他免疫細(xì)胞的表面.PD-1與其在腫瘤細(xì)胞上高度表達(dá)的配體PD-L1相互作用后,發(fā)生一定的構(gòu)象變化,介導(dǎo)胞內(nèi)信號通路從而抑制T細(xì)胞的增殖、活化和細(xì)胞殺傷性功能[7-12].前人的研究已表明,PD-1或PD-L1的基因敲除或抗體抑制可以增強小鼠免疫系統(tǒng)的抗腫瘤功能,這表明阻斷PD-1和PDL1之間的相互作用可能為腫瘤免疫治療提供一種有效的策略[13,14].

PD-1/PD-L1抗體抑制劑藥物的研發(fā)目前取得非常顯著的進(jìn)展.2014年,美國食品藥品監(jiān)督管理局批準(zhǔn)了第一款PD-1抗體Pembrolizumab用于治療晚期黑色素瘤之后,一系列PD-1/PD-L1單克隆抗體被應(yīng)用在了非小細(xì)胞肺癌、肝癌和食管胃交界癌等多種腫瘤疾病的臨床治療中[15-17].然而,隨著應(yīng)用的不斷深入,抗體半衰期長、慢性免疫毒性、組織滲透有限、存儲和運輸成本高等難以避免的缺點逐漸暴露出來[18-20].另一方面,抗體的同質(zhì)化研究過多,造成了較大的資源浪費.小分子抑制劑具備較好的腫瘤組織滲透性、相對穩(wěn)定的生物安全性和良好的口服利用性等優(yōu)勢,已成為PD-1/PD-L1抑制劑藥物的下一個研發(fā)熱點[18,21-23].

百時美施貴寶(Bristol-Myers Squibb,BMS)于2015年公開了一系列非肽基聯(lián)苯類小分子抑制劑,對PD-1/PD-L1結(jié)合具有強大的阻斷抑制活性[24].Holak團(tuán)隊[25]曾報道,BMS化合物通過與PD-L1結(jié)合誘導(dǎo)其發(fā)生二聚化,以間接的方式抑制PD-1/PD-L1相互作用.基于這一機制,其他的公司和學(xué)術(shù)團(tuán)隊開發(fā)了一系列不同骨架的衍生物,如Incyte公司的INCB086550 (NCT04629339/NCT 03762447)[26]、紅日藥業(yè)的IMMH-010 (NCT0434 3859)[27,28]、再極藥業(yè)的MAX-10181 (NCT04122 339)、貝達(dá)藥業(yè)的BPI-371153 (NCT05341557)、歌禮藥業(yè)的ACS61 (NCT05287399)以及和譽生物醫(yī)藥的ABSK043 (NCT04964375)等化合物目前已進(jìn)入到了臨床試驗階段.由于目前PD-L1小分子抑制劑的市場空白,尋找更多樣骨架且具備良好用藥性質(zhì)的化合物仍具備重要意義.

數(shù)據(jù)驅(qū)動的計算方法已成為藥物開發(fā)的重要工具,PD-1/PD-L1小分子抑制劑亦不例外[29].基于結(jié)構(gòu)的方法,如藥效團(tuán)分析、分子對接和分子動力學(xué)(molecular dynamics,MD)模擬,在先前的研究中被廣泛應(yīng)用于靶向PD-L1二聚體的小分子抑制劑的虛擬篩選[29-33].本文基于各種機器學(xué)習(xí)算法和分子描述符或指紋構(gòu)建了一系列基于配體方法的分類和回歸模型,以預(yù)測ZINC15[34]中類藥物化合物對PD-1/PD-L1相互作用的抑制活性.具有高預(yù)測活性的化合物將繼續(xù)通過藥物相似性、藥代動力學(xué)篩選并以分子對接和MD模擬進(jìn)行基于結(jié)構(gòu)方法上的驗證,以最終獲得具備PD-L1抑制潛力的小分子化合物.

2 研究方法

2.1 數(shù)據(jù)獲取及整理

本文用于訓(xùn)練及測試的數(shù)據(jù)集來源于與PDL1小分子抑制劑相關(guān)的37篇研究性論文及16項專利(見補充材料表S1 (online)),為避免因?qū)嶒灱夹g(shù)手段不同而導(dǎo)致化合物對PD-L1抑制活性測定數(shù)據(jù)水平不一致的情況,本文僅收錄了以半抑制濃度(half-maximal inhibitor concentration,IC50)為指標(biāo)的均相時間分辨熒光(homogeneous timeresolved fluorescence,HTRF)的實驗數(shù)據(jù).

參考研究性論文及專利中的結(jié)構(gòu)示意圖,本文利用ChemDraw獲取并檢查化合物的簡化分子線性輸入規(guī)范(simplified molecular input line entry system,SMILES)字符串.IC50不高于1 μmol/L的化合物定義為陽性樣本(即對PD-L1具備抑制活性),而IC50高于10 μmol/L的化合物定義為陰性樣本(即對PD-L1不具備抑制活性),為減小用于分類模型構(gòu)建的兩類樣本數(shù)量不平衡問題,一項細(xì)胞水平的高通量篩選(high throughput screening,HTS)實驗記錄(PubChem BioAssay AID: 2316)部分?jǐn)?shù)據(jù)被引入補充陰性樣本.最后,僅具備明確IC50測定值而非測定范圍的化合物被收錄于回歸模型的數(shù)據(jù)集中,IC50的對數(shù)轉(zhuǎn)化值pIC50(即-lgIC50)作為樣本標(biāo)簽.

2.2 數(shù)據(jù)集聚類

為了確保機器學(xué)習(xí)模型盡可能均勻地獲得數(shù)據(jù)集中不同結(jié)構(gòu)的信息,本文對數(shù)據(jù)集中的化合物先進(jìn)行了聚類處理.首先,2130個分類模型陽性樣本和1099個回歸模型樣本分別轉(zhuǎn)化為600和350種僅保留環(huán)形結(jié)構(gòu)和連接環(huán)形結(jié)構(gòu)的最短路徑的Murcko骨架[35].分子骨架以2為最大半徑,2048為向量長度轉(zhuǎn)化為圓形擴(kuò)展指紋(extendedconnectivity finger print,ECFP)[36]后,以平均為鏈接算法、歐氏距離為計算方式對兩類骨架的ECFP進(jìn)行分層聚類,簇的數(shù)量由5—20之間對應(yīng)的最佳輪廓系數(shù)決定,后續(xù)模型訓(xùn)練將從簇內(nèi)分層抽樣進(jìn)行數(shù)據(jù)集劃分(圖1).

圖1 數(shù)據(jù)集聚類 (a)分類模型陽性樣本骨架的14個簇;(b)回歸模型樣本骨架的13個簇.環(huán)形樹狀圖的顏色代表不同的分子結(jié)構(gòu)骨架簇,顏色越接近,骨架越相似,一旁的化學(xué)結(jié)構(gòu)圖是每個簇中最具代表性的骨架Fig.1.Dataset clustering: (a) 14 clusters of scaffolds of active compounds in the classification models;(b) 13 clusters of scaffolds of compounds in the regression models.The colors of the circular dendrograms represent different clusters of scaffolds,and the closer the colors are,the more similar the scaffolds are.The chemical structure diagrams on the side are the most representative scaffolds of each cluster.

2.3 分子表征

由ChemDraw得到的化合物SMILES字符串分別轉(zhuǎn)化為三類分子描述符(RDKit,PaDEL 1D&2D,Mordred)和三類分子指紋(ECFP,MACCS,PubChem)以作為化合物的特征向量.開源Java軟件PaDEL-Descriptor v2.0用于計算PaDEL 1D&2D描述符[37]、MACCS分子指紋[38]和Pub Chem分子指紋[37].RDKit描述符和ECFP分子指紋[36]均由化學(xué)信息處理程序包RDKit轉(zhuǎn)化,其中ECFP分子指紋為2048維向量,指示最大以兩個原子為半徑的結(jié)構(gòu)碎片存在與否.此外,描述符計算軟件Mordred[39]被用于最后一類特征向量的轉(zhuǎn)化.

2.4 機器學(xué)習(xí)

分別使用五種機器學(xué)習(xí)算法構(gòu)建PD-L1小分子抑制劑分類模型和抑制活性IC50預(yù)測回歸模型,其中邏輯回歸(logistic regression,LR)、鄰近算法(K-nearest neighbor,KNN)、支持向量機(support vector machine,SVM)、隨機森林(random forest,RF)和多層感知機(multilayer perceptron,MLP)用于分類任務(wù)模型構(gòu)建,SVM、嶺回歸(ridge regression)、高斯過程回歸(Gaussian process regression,GPR)、RF和MLP用于回歸任務(wù)模型構(gòu)建.兩類模型數(shù)據(jù)集的80%為訓(xùn)練數(shù)據(jù)集,另外20%為測試數(shù)據(jù)集.所有分子描述符及分子指紋的特征值刪除空缺之后,基于訓(xùn)練集進(jìn)行min-max標(biāo)準(zhǔn)化處理.各類算法的最佳超參數(shù)由基于Scikitlearn網(wǎng)格搜索方法的五重交叉驗證(5-fold cross validation)確定,此過程中,分類模型以馬修斯相關(guān)系數(shù)(matthews correlation coefficient,MCC),回歸模型以預(yù)測值的平均絕對誤差(mean absolute error,MAE)為超參數(shù)選擇標(biāo)準(zhǔn).

2.5 模型評價標(biāo)準(zhǔn)

分類模型性能由靈敏度(sensitivity,SE)、特異度(specificity,SP)、準(zhǔn)確度(accuracy,ACC)、馬修斯相關(guān)系數(shù)(Matthews correlation coefficient,MCC)進(jìn)行泛化評估,其計算公式如下:

其中,TP(true positive)為分類模型判定為具備抑制活性且實驗結(jié)果亦為具備抑制活性的化合物樣本數(shù),FP(false positive)為分類模型判定為具備抑制活性但實驗結(jié)果為不具備抑制活性的化合物樣本數(shù),TN(true negative)為分類模型判定為不具備抑制活性且實驗結(jié)果亦為不具備抑制活性的化合物樣本數(shù),FN(false negative)為分類模型判定為不具備抑制活性但實驗結(jié)果為具備抑制活性的化合物樣本數(shù).

回歸模型性能由平均絕對誤差(mean absolute error,MAE)、均方根誤差(root mean-square error,RMSE)和決定系數(shù)(correlation coefficient,R2)進(jìn)行泛化評估,其計算公式如下:

其中,n為回歸模型待評估數(shù)據(jù)集樣本數(shù),Ytrue為化合物pIC50實驗測定值,Ypred為回歸模型對化合物pIC50的預(yù)測估計值.

2.6 藥物相似性及ADMET檢驗

本文選取來自ZINC15數(shù)據(jù)庫[34]的7400926個可商業(yè)購買(截至2023年2月24日)、電中性且收錄三維結(jié)構(gòu)信息的類藥小分子作為候選化合物篩選池,盡管這些分子已通過ZINC15的藥物相似性檢驗,但由于本文所篩選的化合物針對蛋白質(zhì)相互作用,傳統(tǒng)的藥物相似性指標(biāo)(quantitative estimate of drug-likeness,QED)已經(jīng)不再適用[40].Kosugi和Ohue[41]提出了一種更加適合于針對蛋白質(zhì)相互作用抑制劑篩選的藥物相似性指標(biāo)(quantitative estimate index for compounds targeting protein-protein interactions,QEPPI),化合物的QEPPI分?jǐn)?shù)為一個介于0—1之間的值,數(shù)值越大代表化合物的蛋白質(zhì)相互作用抑制劑藥物相似性越高,本文以0.7為閾值篩選高類藥性化合物.

ADMET代表化合物的吸收(absorption)、分配(distribution)、代謝(metabolism)、排泄(excretion)和毒性(toxicity)等重要用藥性質(zhì),在早期藥物篩選中十分重要,本文采用ADMETlab 2.0對化合物進(jìn)行ADMET篩選[42].

2.7 分子對接

基于結(jié)構(gòu)的活性預(yù)測方法需要蛋白質(zhì)靶點的晶體結(jié)構(gòu)信息,本文選取的PD-L1二聚體結(jié)構(gòu)來源于蛋白質(zhì)數(shù)據(jù)庫PDB(ID: 7DY7).本文選擇PD-L1二聚體中的結(jié)合界面區(qū)域為對接結(jié)合口袋,基于如下理由: 1)前期的研究表明PD-L1二聚體中的結(jié)合界面區(qū)域是很有潛力的藥物結(jié)合位點[43,44];2)分子對接結(jié)果顯示小分子在該區(qū)域的結(jié)合親和力最強.

分子對接利用AutoDock Vina進(jìn)行[45],Auto-DockTools將PD-L1二聚體轉(zhuǎn)化為pdbqt格式,對接網(wǎng)格盒子尺寸為40 ?×30 ?×40 ?,中心坐標(biāo)為(144.799 ?,-9.364 ?,16.163 ?).每個化合物各隨機生成50種對接姿態(tài)和位置,根據(jù)對接得分進(jìn)行排名,最佳得分前十位小分子對應(yīng)的結(jié)合構(gòu)象將用于分子動力學(xué)模擬的初始構(gòu)象.此外,為驗證通過上述流程得到的候選化合物作用機制是否與前人的研究結(jié)果相近,本文利用LigPlot+研究小分子與PD-L1二聚體的相互作用[46].

2.8 分子動力學(xué)模擬

本文所有的分子動力學(xué)模擬采用CHARMM 36全原子力場[47,48],化合物配體的力場參數(shù)通過CGenFF程序獲取[47-50].十個高對接分?jǐn)?shù)的候選化合物和兩個對照化合物(BMS202以及INCB086550)與PD-L1二聚體復(fù)合系統(tǒng)搭建于10nm×10nm×10nm盒子內(nèi),每個盒子填充水分子并以NaCl中和體系電荷數(shù),整個搭建過程由GROMACS[51]工具gmx solvate和gmx genion實現(xiàn).所有的分子動力學(xué)模擬工作均使用GROMACS 2019.6程序包運行,模擬步長為2.0 fs,采用等溫等壓(NPT)系綜,模擬時長為100 ns.模擬盒在x軸、y軸和z軸3個方向上均設(shè)定了周期性邊界條件.采用半各向同性的Parrinello-Rahman[52]方法將系統(tǒng)壓強維持在1 bar (1 bar=105Pa),壓縮系數(shù)設(shè)定為4.5×10-5,弛豫時間為5 ps.此外,系統(tǒng)采用Nose-Hoover[53,54]控溫方法將配體-蛋白質(zhì)復(fù)合物和溶劑進(jìn)行耦合,溫度維持在310 K,弛豫時間為1 ps.長距離靜電相互作用使用Particle-Mesh Ewald (PME)[55]方法計算,短距離靜電相互作用的截止距離設(shè)置為1.2 nm,LINCS (LINear constraint solver)[56]算法用于約束含H原子的鍵長.

2.9 MM/PBSA結(jié)合自由能計算

本文以gmx_MMPBSA[57]為工具利用分子力學(xué)泊松-玻爾茲曼表面積法(molecular mechanics Poisson-Boltzmann surface area,MM/PBSA)來計算各個體系最后10 ns軌跡的平均結(jié)合自由能(ΔG),其計算公式如下:

其中Gcomplex為蛋白質(zhì)-配體復(fù)合物的自由能,Gprotein和Gligand分別為蛋白質(zhì)和配體各自的自由能.各組分的自由能計算公式如下:

其中EMM為氣相結(jié)合能量,由范德瓦耳斯項Evdw和靜電項Eele構(gòu)成;Gsol為溶劑化自由能,由極性EPB和非極性ESA貢獻(xiàn)構(gòu)成;TΔS為熵變,因其計算成本較高且未必能夠提高自由能的精度,本文暫不對其進(jìn)行計算.

3 研究結(jié)果

3.1 分類模型性能

30種分類模型分別基于LR,KNN,SVM,RF和MLP五種算法以及RDKit,PaDEL 1D&2D,Mordred分子描述符和ECFP,PubChem,MACCS分子指紋6種分子表征方式所建立.圖2為各個模型經(jīng)五折驗證并設(shè)置最佳超參數(shù)后在測試集上的性能表現(xiàn),其中以ECFP為輸入的KNN模型表現(xiàn)出最佳性能,SE(陽性樣本預(yù)測正確率)=0.9937,SP(陰性樣本預(yù)測正確率)=0.9781,ACC(全部樣本預(yù)測正確率)=0.9859,MCC(綜合評價兩類樣本預(yù)測性能指標(biāo))=0.9720.考慮到分子描述符或分子指紋的計算較為耗時,同時獲取740萬余個類藥小分子的全部六種輸入向量的計算量更大,因此,僅采用以ECFP為輸入用于后續(xù)的分類篩選.另一方面,相對于LR和MLP兩種模型,KNN,SVM和RF模型對于可能具備噪聲數(shù)據(jù)的任務(wù)具有較好的魯棒性,因此,選擇KNN,SVM和RF模型作為分類篩選的模型;當(dāng)有兩種或以上的模型支持小分子為活性時,則認(rèn)定該化合物對PD-L1具備抑制活性.

圖2 分類模型性能表現(xiàn) (a)靈敏度;(b)特異度;(c)準(zhǔn)確度;(d)馬修斯相關(guān)系數(shù)Fig.2.Performance of binary models for classification tasks: (a) SE;(b) SP;(c) ACC;(d) MCC.

3.2 分類模型解釋

為了進(jìn)一步探索輸入向量的特征與分類之間的相關(guān)性,本文用活性和非活性化合物之間的RDKit,PaDEL和Mordred三種描述符來分析基于特征重要性而構(gòu)建的RF模型中權(quán)重最高的5個特征的分布差異(見補充材料圖S1 (online)).盡管部分特征在兩類樣本間的分布差異是十分顯著的,但由于分子描述符的特征信息難以解釋,特征值根據(jù)理化性質(zhì)及矩陣運算得出,其大小本身不具備具象含義,我們的認(rèn)知只能停留在較淺薄的層面.

與描述符不同的是,分子指紋ECFP的位點指示局部亞結(jié)構(gòu)的存在或不存在,這可能揭示分子結(jié)構(gòu)與對PD-1/PD-L1的抑制活性之間的關(guān)系.分別統(tǒng)計了前十位活性化合物中比例顯著高于非活性化合物的子結(jié)構(gòu)和前十位活性化合物中比例顯著低于非活性化合物的子結(jié)構(gòu)(圖3).帶有黃色芳香性原子的環(huán)狀結(jié)構(gòu)在活性化合物中的計數(shù)遠(yuǎn)多于其在非活性化合物中的計數(shù),這意味著芳香性結(jié)構(gòu)片段可能對化合物的PD-L1抑制活性有正向貢獻(xiàn),相反,在非活性化合物中常出現(xiàn)的帶灰色原子的脂肪鏈片段則對活性沒有正向貢獻(xiàn).

圖3 分子結(jié)構(gòu)片段在兩類樣本間的計數(shù)差異 (a)活性化合物計數(shù)占優(yōu)的結(jié)構(gòu)片段;(b)非活性化合物計數(shù)占優(yōu)的結(jié)構(gòu)片段.結(jié)構(gòu)片段的中心原子以藍(lán)色突出顯示;芳香性原子被著色為黃色,而脂肪烴鏈原子則被著色為灰色Fig.3.Count difference of fragments of structures between active and inactive compounds: (a) Fragments in active compounds with a proportion higher than inactive compounds;(b) fragments in active compounds with a proportion lower than inactive compounds.The center atoms of substructure fragments are highlighted in blue;aromaticity atoms are colored in yellow and aliphatic hydrocarbon atoms are colored in gray.

值得注意的是,ECFP985和ECFP1161片段清楚地表示了化合物的聯(lián)苯特性,這也是BMS化合物的核心結(jié)構(gòu)特征之一[24,58].此外,苯甲基(ECFP253)、吡啶(ECFP1453)和與苯環(huán)相連的醚鍵(ECFP1971)通常存在于多數(shù)表現(xiàn)出PD-L1抑制活性的小分子結(jié)構(gòu)中.前人的分子對接分析表明,苯甲基能夠與PD-L1的Ala121和Met115產(chǎn)生強烈的疏水相互作用[59].Lu等[60]發(fā)現(xiàn)PD-L1-BMS202復(fù)合物中,吡啶環(huán)的氮原子周圍有很大的空間,因此他們在該位點附加了一系列的取代基,以提高配體的結(jié)合親和力.由于該片段的重要性,許多其他研究團(tuán)隊也將工作重點放在吡啶的修飾上[59,61,62].除了環(huán)狀結(jié)構(gòu),連接環(huán)狀結(jié)構(gòu)的路徑可能也是值得關(guān)注的組成部分,其中,以醚鍵連接六元環(huán)或五元環(huán)的活性化合物約占50%[43,63-65].

3.3 回歸模型性能

30種回歸模型分別基于SVM,Ridge Regression,GPR,RF和MLP五種算法以及RDKit,PaDEL 1D&2D,Mordred分子描述符和ECFP,PubChem,MACCS分子指紋6種分子表征方式所建立.圖4為各個模型經(jīng)五折驗證后、最佳超參數(shù)設(shè)置下在測試集上的性能表現(xiàn),其中以ECFP為輸入的SVM模型以MAE=0.4503,RMSE=0.6375,GPR模型以MAE=0.4557,RMSE=0.6375的性能表現(xiàn)顯著優(yōu)于其他模型.圖5展示了兩種模型在訓(xùn)練集及測試集所有樣本點的預(yù)測結(jié)果及偏差,絕大部分的樣本點預(yù)測偏差在±1范圍內(nèi),以ECFP為輸入的SVM模型在測試上的決定系數(shù)R2=0.782,以ECFP為輸入的GPR模型在測試上的決定系數(shù)R2=0.787,兩類模型預(yù)測值的平均值將作為候選化合物的PD-L1抑制活性預(yù)測值.

圖4 回歸模型性能表現(xiàn) (a)平均絕對誤差;(b)均方根誤差Fig.4.Performance of continuous models for regression tasks: (a) MAE;(b) RMSE.

圖5 兩種最佳性能回歸模型預(yù)測結(jié)果 (a),(c)樣本預(yù)測值與標(biāo)簽值分布;(b),(d)樣本預(yù)測偏差Fig.5.Prediction results of two best-performing continuous models: (a),(c) Distribution of predicted values and label values;(b),(d) prediction residuals.

3.4 虛擬篩選

ZINC15數(shù)據(jù)庫中電中性且收錄三維結(jié)構(gòu)信息的類藥物小分子構(gòu)成本文的化合物篩選池(7400926個小分子).整個虛擬篩選過程包含4個環(huán)節(jié): 1)同時使用以ECFP為輸入的KNN,SVM和RF等3種分類模型,其中至少有2種判定結(jié)果為對PD-L1蛋白具有抑制活性,則認(rèn)定該分子具有抑制活性.2)同時使用以ECFP為輸入的SVM和GPR回歸模型預(yù)測pIC50值,兩模型pIC50的平均值小于7 (IC50 ＞ 100 nmol/L)的化合物將被剔除.3)計算小分子化合物的QEPPI分?jǐn)?shù)(見補充材料表S2 (online)),并保留QEPPI得分超過0.7的候選化合物.盡管初始篩選池的化合物在ZINC15中根據(jù)Lipinski五原則[66]被定義為類藥化合物,但其中部分化合物因分子量過小可能難以充分結(jié)合到PD-L1二聚體的相互作用界面,因此,重新評估他們的藥物相似性仍具有重要意義.4)使用ADMETlab 2.0[42]計算小分子的AD MET性質(zhì),以進(jìn)一步篩選具有應(yīng)用潛能的小分子(見補充材料表S3 (online)).最終,通過整個虛擬篩選流程,42個候選化合物被認(rèn)為對PD-L1具有較強的抑制活性并具備良好的藥用性能(圖6).

圖6 虛擬篩選流程Fig.6.Workflow of virtual screening.

3.5 分子對接

為了驗證篩選結(jié)果并通過基于結(jié)構(gòu)的方法進(jìn)一步識別有抑制潛力的化合物,利用AutoDock Vina將42個候選化合物以及BMS202和INCB-086550對接到PD-L1二聚體的IgV結(jié)構(gòu)域(PDB ID:7DY7).具有最低對接評分的10種化合物被認(rèn)為有潛力的PD-L1抑制劑并繼續(xù)用于之后的分析(表1),其中值得注意的是,10種候選化合物對接分?jǐn)?shù)均低于臨床II期試驗化合物INCB086550的對接分?jǐn)?shù),因此可以認(rèn)為先前構(gòu)建的分類和回歸模型是篩選PD-L1小分子抑制劑的有效工具.分析12種化合物與PD-L1二聚體的相互作用可知(見補充材料圖S2 (online)),TYR56(A),TYR56(B),MET115(A),MET115(B),ALA121(A)和TYR123(A)(括弧中A或B分別表示PD-L1二聚體中的單體A或B)等氨基酸殘基與所有配體均發(fā)生了較強的疏水相互作用,這意味著這些小分子的結(jié)合位置和結(jié)合模式可能是近似的.此外,除了ZINC000021874692和ZINC000021874694這一對化合物為旋光異構(gòu)體,其余化合物的結(jié)構(gòu)差異較大,相似度較低,意味著其可為未來的濕實驗篩選提供較高的容錯空間(見補充材料圖S3 (online)).

表1 最低對接評分的10種候選化合物及BMS202和INCB086550的對接結(jié)果Table 1.Docking results for the top 10 hits with BMS202 and INCB086550.

3.6 分子動力學(xué)模擬

為驗證配體與PD-L1二聚體結(jié)合的穩(wěn)定性,計算了12種配體-蛋白質(zhì)復(fù)合物100 ns軌跡中蛋白質(zhì)骨架和配體小分子構(gòu)象基于蛋白質(zhì)骨架疊合的均方根偏差(root mean square deviation,RMSD)(圖7).所有體系的PD-L1二聚體骨架和配體的RMSD在50 ns后均穩(wěn)定在0.8 nm以下,模擬軌跡達(dá)到平衡,蛋白質(zhì)與配體結(jié)合穩(wěn)定.值得注意的是,以INCB086550為例的分子量較大的小分子與以ZINC000019770413為例的分子量較小的小分子相比,配體構(gòu)象RMSD波動較大.由圖8的蛋白質(zhì)-配體結(jié)合模式可知,起到關(guān)鍵作用的結(jié)構(gòu)片段集中在配體的一端,而另一端暴露于結(jié)合口袋之外.若配體的分子量較大,分子骨架較長,其在結(jié)合區(qū)域外的部分更多,這部分運動更為自由,這可能是部分小分子構(gòu)象RMSD波動相對較大的原因.

圖7 RMSD (a) PD-L1二聚體骨架;(b)配體Fig.7.RMSD: (a) Backbone of PD-L1 dimer;(b) ligand.

圖8 配體與PD-L1結(jié)合模式和關(guān)鍵殘基 (a) ZINC000021 723762;(b) ZINC000021 874692;(c) ZINC000175 468610;(d) ZINC 000019 770413;(e) ZINC000021 874694;(f) ZINC000952 973550;(g) ZINC000064 987401;(h) ZINC000003 908573;(i) ZINC 000020 538424;(j) ZINC000004 063088;(k) BMS202;(l) INCB086550Fig.8.Ligand binding mode with PD-L1 and key residues: (a) ZINC000021 723762;(b) ZINC000021 874692;(c) ZINC000175 468610;(d) ZINC000019 770413;(e) ZINC000021 874694;(f) ZINC000952 973550;(g) ZINC000064 987401;(h) ZINC000003 908573;(i) ZINC000020 538424;(j) ZINC000004 063088;(k) BMS202;(l) NCB086550.

3.7 MM/PBSA結(jié)合自由能計算

分子動力學(xué)模擬的最后10 ns軌跡被用于MM/PBSA計算以評估配體小分子與PD-L1二聚體間的結(jié)合作用強度.由表2可知,大部分候選化合物的結(jié)合自由能均與對照組的結(jié)合自由能處于同一水平,而ZINC000021723762,ZINC0000218 74694和ZINC000004063088甚至顯著優(yōu)于對照組水平,范德瓦耳斯相互作用為驅(qū)動這三類化合物區(qū)別于其他化合物的主要因素.

表2 結(jié)合自由能計算結(jié)果Table 2.Results of MMPBSA.

將結(jié)合自由能分解至每個氨基酸殘基與配體分子的相互作用上,提取了每個體系中對結(jié)合強度貢獻(xiàn)度最大的5個關(guān)鍵殘基.如圖8所示,TYR56(B),TYR123(A),MET115(B),ALA121(A)在絕大部分體系中都起到關(guān)鍵作用,這意味著12種配體結(jié)合于PD-L1二聚體相互作用交界面的模式是近似的.前人研究揭示出,MET115(B)和ALA 121(A)與化合物中的芳香環(huán)能夠發(fā)生強烈疏水相互作用,而TYR56(B)能夠與化合物的苯環(huán)形成π-π堆積以增強結(jié)合穩(wěn)定性[32,33,59,67],這些作用在結(jié)合自由能貢獻(xiàn)中均得到了體現(xiàn).

4 結(jié)論

通過從相關(guān)研究文獻(xiàn)及專利收集PD-L1小分子抑制活性數(shù)據(jù)集,并對數(shù)據(jù)集內(nèi)化合物以結(jié)構(gòu)相似性進(jìn)行分層抽樣,本文構(gòu)建了各30種基于不同分子表征方法和算法的機器學(xué)習(xí)活性判定分類模型和活性強度回歸預(yù)測模型,其中性能最佳分類模型分類正確率可達(dá)98%以上,回歸模型預(yù)測平均絕對誤差在0.5以下,具備良好的應(yīng)用價值.將以上模型并結(jié)合藥用性質(zhì)篩選工具構(gòu)成的虛擬篩選方法應(yīng)用于ZINC15大型小分子數(shù)據(jù)庫,獲得了42種有潛力的PD-L1小分子抑制劑,其中的10種通過分子對接、分子動力學(xué)模擬和結(jié)合自由能估計驗證發(fā)現(xiàn),候選化合物的作用機制與對照化合物相近,部分化合物的結(jié)合強度甚至優(yōu)于對照化合物,這意味著本文前期建立的機器學(xué)習(xí)模型是可以幫助加速PD-L1小分子抑制劑虛擬篩選的有效工具.

然而,完整的藥物研發(fā)并不是一個簡單的工程,本文僅使用計算方法幫助加速PD-L1小分子抑制劑的研發(fā),在完整的研發(fā)產(chǎn)業(yè)鏈中處于較上游的位置.為盡可能使得本文的篩選結(jié)果得到更進(jìn)一步的有效認(rèn)證,后續(xù)分子水平、細(xì)胞水平以及動物模型水平的濕實驗驗證仍是必不可缺的,計算機技術(shù)目前僅能扮演錦上添花的角色.相信隨著生物計算領(lǐng)域與生物技術(shù)領(lǐng)域的合作加深,藥物研發(fā)將能夠向一個盡可能高效、低成本的方向發(fā)展.