靶向PD-L1蛋白的計(jì)算機(jī)輔助藥物篩選*

林開(kāi)東 林曉倩2) 林緒波?

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

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

1 引言

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

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

百時(shí)美施貴寶(Bristol-Myers Squibb,BMS)于2015年公開(kāi)了一系列非肽基聯(lián)苯類小分子抑制劑,對(duì)PD-1/PD-L1結(jié)合具有強(qiáng)大的阻斷抑制活性[24].Holak團(tuán)隊(duì)[25]曾報(bào)道,BMS化合物通過(guò)與PD-L1結(jié)合誘導(dǎo)其發(fā)生二聚化,以間接的方式抑制PD-1/PD-L1相互作用.基于這一機(jī)制,其他的公司和學(xué)術(shù)團(tuán)隊(duì)開(kāi)發(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ù)生物醫(yī)藥的ABSK043 (NCT04964375)等化合物目前已進(jìn)入到了臨床試驗(yàn)階段.由于目前PD-L1小分子抑制劑的市場(chǎng)空白,尋找更多樣骨架且具備良好用藥性質(zhì)的化合物仍具備重要意義.

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

2 研究方法

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

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

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

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

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

圖1 數(shù)據(jù)集聚類 (a)分類模型陽(yáng)性樣本骨架的14個(gè)簇;(b)回歸模型樣本骨架的13個(gè)簇.環(huán)形樹(shù)狀圖的顏色代表不同的分子結(jié)構(gòu)骨架簇,顏色越接近,骨架越相似,一旁的化學(xué)結(jié)構(gòu)圖是每個(gè)簇中最具代表性的骨架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)以作為化合物的特征向量.開(kāi)源Java軟件PaDEL-Descriptor v2.0用于計(jì)算PaDEL 1D&2D描述符[37]、MACCS分子指紋[38]和Pub Chem分子指紋[37].RDKit描述符和ECFP分子指紋[36]均由化學(xué)信息處理程序包RDKit轉(zhuǎn)化,其中ECFP分子指紋為2048維向量,指示最大以兩個(gè)原子為半徑的結(jié)構(gòu)碎片存在與否.此外,描述符計(jì)算軟件Mordred[39]被用于最后一類特征向量的轉(zhuǎn)化.

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

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

2.5 模型評(píng)價(jià)標(biāo)準(zhǔn)

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

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

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

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

2.6 藥物相似性及ADMET檢驗(yàn)

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

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

2.7 分子對(duì)接

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

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

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

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

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

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

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

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

3 研究結(jié)果

3.1 分類模型性能

30種分類模型分別基于LR,KNN,SVM,RF和MLP五種算法以及RDKit,PaDEL 1D&2D,Mordred分子描述符和ECFP,PubChem,MACCS分子指紋6種分子表征方式所建立.圖2為各個(gè)模型經(jīng)五折驗(yàn)證并設(shè)置最佳超參數(shù)后在測(cè)試集上的性能表現(xiàn),其中以ECFP為輸入的KNN模型表現(xiàn)出最佳性能,SE(陽(yáng)性樣本預(yù)測(cè)正確率)=0.9937,SP(陰性樣本預(yù)測(cè)正確率)=0.9781,ACC(全部樣本預(yù)測(cè)正確率)=0.9859,MCC(綜合評(píng)價(jià)兩類樣本預(yù)測(cè)性能指標(biāo))=0.9720.考慮到分子描述符或分子指紋的計(jì)算較為耗時(shí),同時(shí)獲取740萬(wàn)余個(gè)類藥小分子的全部六種輸入向量的計(jì)算量更大,因此,僅采用以ECFP為輸入用于后續(xù)的分類篩選.另一方面,相對(duì)于LR和MLP兩種模型,KNN,SVM和RF模型對(duì)于可能具備噪聲數(shù)據(jù)的任務(wù)具有較好的魯棒性,因此,選擇KNN,SVM和RF模型作為分類篩選的模型;當(dāng)有兩種或以上的模型支持小分子為活性時(shí),則認(rèn)定該化合物對(duì)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三種描述符來(lái)分析基于特征重要性而構(gòu)建的RF模型中權(quán)重最高的5個(gè)特征的分布差異(見(jiàn)補(bǔ)充材料圖S1 (online)).盡管部分特征在兩類樣本間的分布差異是十分顯著的,但由于分子描述符的特征信息難以解釋,特征值根據(jù)理化性質(zhì)及矩陣運(yùn)算得出,其大小本身不具備具象含義,我們的認(rèn)知只能停留在較淺薄的層面.

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

圖3 分子結(jié)構(gòu)片段在兩類樣本間的計(jì)數(shù)差異 (a)活性化合物計(jì)數(shù)占優(yōu)的結(jié)構(gòu)片段;(b)非活性化合物計(jì)數(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)中.前人的分子對(duì)接分析表明,苯甲基能夠與PD-L1的Ala121和Met115產(chǎn)生強(qiáng)烈的疏水相互作用[59].Lu等[60]發(fā)現(xiàn)PD-L1-BMS202復(fù)合物中,吡啶環(huán)的氮原子周圍有很大的空間,因此他們?cè)谠撐稽c(diǎn)附加了一系列的取代基,以提高配體的結(jié)合親和力.由于該片段的重要性,許多其他研究團(tuán)隊(duì)也將工作重點(diǎ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為各個(gè)模型經(jīng)五折驗(yàn)證后、最佳超參數(shù)設(shè)置下在測(cè)試集上的性能表現(xiàn),其中以ECFP為輸入的SVM模型以MAE=0.4503,RMSE=0.6375,GPR模型以MAE=0.4557,RMSE=0.6375的性能表現(xiàn)顯著優(yōu)于其他模型.圖5展示了兩種模型在訓(xùn)練集及測(cè)試集所有樣本點(diǎn)的預(yù)測(cè)結(jié)果及偏差,絕大部分的樣本點(diǎn)預(yù)測(cè)偏差在±1范圍內(nèi),以ECFP為輸入的SVM模型在測(cè)試上的決定系數(shù)R2=0.782,以ECFP為輸入的GPR模型在測(cè)試上的決定系數(shù)R2=0.787,兩類模型預(yù)測(cè)值的平均值將作為候選化合物的PD-L1抑制活性預(yù)測(cè)值.

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

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

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

3.5 分子對(duì)接

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

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

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

為驗(yàn)證配體與PD-L1二聚體結(jié)合的穩(wěn)定性,計(jì)算了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波動(dòng)較大.由圖8的蛋白質(zhì)-配體結(jié)合模式可知,起到關(guān)鍵作用的結(jié)構(gòu)片段集中在配體的一端,而另一端暴露于結(jié)合口袋之外.若配體的分子量較大,分子骨架較長(zhǎng),其在結(jié)合區(qū)域外的部分更多,這部分運(yùn)動(dòng)更為自由,這可能是部分小分子構(gòu)象RMSD波動(dòng)相對(duì)較大的原因.

圖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é)合自由能計(jì)算

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

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

將結(jié)合自由能分解至每個(gè)氨基酸殘基與配體分子的相互作用上,提取了每個(gè)體系中對(duì)結(jié)合強(qiáng)度貢獻(xiàn)度最大的5個(gè)關(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ā)生強(qiáng)烈疏水相互作用,而TYR56(B)能夠與化合物的苯環(huán)形成π-π堆積以增強(qiáng)結(jié)合穩(wěn)定性[32,33,59,67],這些作用在結(jié)合自由能貢獻(xiàn)中均得到了體現(xiàn).

4 結(jié)論

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

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