靶向PD-L1蛋白的計算機輔助藥物篩選*

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

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

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

1 引言

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

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

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

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

2 研究方法

2.1 數據獲取及整理

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

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

2.2 數據集聚類

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

圖1 數據集聚類 (a)分類模型陽性樣本骨架的14個簇;(b)回歸模型樣本骨架的13個簇.環形樹狀圖的顏色代表不同的分子結構骨架簇,顏色越接近,骨架越相似,一旁的化學結構圖是每個簇中最具代表性的骨架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字符串分別轉化為三類分子描述符(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]均由化學信息處理程序包RDKit轉化,其中ECFP分子指紋為2048維向量,指示最大以兩個原子為半徑的結構碎片存在與否.此外,描述符計算軟件Mordred[39]被用于最后一類特征向量的轉化.

2.4 機器學習

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

2.5 模型評價標準

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

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

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

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

2.6 藥物相似性及ADMET檢驗

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

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

2.7 分子對接

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

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

2.8 分子動力學模擬

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

2.9 MM/PBSA結合自由能計算

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

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

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

3 研究結果

3.1 分類模型性能

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

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

3.2 分類模型解釋

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

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

圖3 分子結構片段在兩類樣本間的計數差異 (a)活性化合物計數占優的結構片段;(b)非活性化合物計數占優的結構片段.結構片段的中心原子以藍色突出顯示;芳香性原子被著色為黃色,而脂肪烴鏈原子則被著色為灰色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片段清楚地表示了化合物的聯苯特性,這也是BMS化合物的核心結構特征之一[24,58].此外,苯甲基(ECFP253)、吡啶(ECFP1453)和與苯環相連的醚鍵(ECFP1971)通常存在于多數表現出PD-L1抑制活性的小分子結構中.前人的分子對接分析表明,苯甲基能夠與PD-L1的Ala121和Met115產生強烈的疏水相互作用[59].Lu等[60]發現PD-L1-BMS202復合物中,吡啶環的氮原子周圍有很大的空間,因此他們在該位點附加了一系列的取代基,以提高配體的結合親和力.由于該片段的重要性,許多其他研究團隊也將工作重點放在吡啶的修飾上[59,61,62].除了環狀結構,連接環狀結構的路徑可能也是值得關注的組成部分,其中,以醚鍵連接六元環或五元環的活性化合物約占50%[43,63-65].

3.3 回歸模型性能

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

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

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

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

3.5 分子對接

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

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

3.6 分子動力學模擬

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

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

圖8 配體與PD-L1結合模式和關鍵殘基 (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結合自由能計算

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

表2 結合自由能計算結果Table 2.Results of MMPBSA.

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

4 結論

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

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