絲氨酸及苯基異絲氨酸類SARS病毒3CL蛋白酶抑制劑HQSAR研究及分子設(shè)計(jì)

何安朕，張繼川，何嚴(yán)萍，王月平**

(1.昆明理工大學(xué) 理學(xué)院應(yīng)用化學(xué)系，云南 昆明 650500；2.云南大學(xué) 教育部自然資源藥物化學(xué)重點(diǎn)實(shí)驗(yàn)室，化學(xué)科學(xué)與工程學(xué)院·藥學(xué)院，云南 昆明 650500)

冠狀病毒是自然界廣泛存在的一種病毒，在系統(tǒng)分類上屬套式病毒目(Nidovirales)冠狀病毒科(Coronaviridae)冠狀病毒屬(Coronavirus)[1].此類病毒是具囊膜(Envelope)、基因組為線性單股正鏈的RNA病毒.目前已知的可感染人類的冠狀病毒一共有7種，分別為HCoV-229E、HCoV-OC43、HCoVNL63、HCoV-HKU1、SARS-CoV、MERS-CoV和SARS-CoV-2. 其中，后3種對(duì)人類健康及公共安全具有嚴(yán)重威脅：SARS-CoV導(dǎo)致了2003年重癥急性呼吸綜合征的爆發(fā)、MERS-CoV引發(fā)了2012年的中東呼吸綜合征，而2019年的新型冠狀病毒肺炎(COVID-19)[2]則因感染SARS-CoV-2所致.由于新冠疫情的全球蔓延，研發(fā)SARS-CoV-2的藥物相關(guān)藥物刻不容緩.鑒于SARS-CoV-2與SARSCoV的基因組具有82%的相似度[3].兩者高度同源，說明這兩種病毒可能均起源于蝙蝠[4].因此，從抗SARS-CoV化合物中發(fā)現(xiàn)SARS-CoV-2抑制劑，是目前研發(fā)治療COVID-19藥物的重要策略之一[5-7].

冠狀病毒基因組編碼的病毒蛋白主要有以下3種：結(jié)構(gòu)蛋白、復(fù)制酶多聚蛋白及附屬蛋白.就SARS冠狀病毒而言，其結(jié)構(gòu)蛋白主要由N蛋白、M蛋白、E蛋白以及S蛋白共同構(gòu)成；附屬蛋白則包括由orfla、orflb等開放閱讀框編碼的各種蛋白質(zhì)[8]；其復(fù)制酶多聚蛋白pp1a和pp1b則由開放閱讀框1a/b編碼，生成后被3CL蛋白酶和類木瓜蛋白酶(PLpro)水解，裁剪成具有獨(dú)立作用的16個(gè)功能性多肽（蛋白）.只有在這些非結(jié)構(gòu)蛋白的參與下，病毒才能進(jìn)行轉(zhuǎn)錄和復(fù)制[9].因此，3CL蛋白酶是抗SARS病毒藥物設(shè)計(jì)的重要靶標(biāo).

3CL蛋白酶的催化活性中心由半胱氨酸(Cys145)和組氨酸(His41)二聯(lián)體構(gòu)成.其中，Cys145為親核進(jìn)攻基團(tuán)，His41則為酸-堿催化殘基[10-11].2003年Rao課題組[12]報(bào)道了肽類抑制劑與3CL蛋白酶結(jié)合的晶體結(jié)構(gòu)，為后續(xù)的藥物設(shè)計(jì)指明了方向.Kenichi Akaji課題組基于3CL蛋白酶晶體結(jié)構(gòu)進(jìn)行合理設(shè)計(jì)，合成了高活性的醛基四肽[7,13].構(gòu)效關(guān)系分析表明，當(dāng)抑制劑能與3CL酶結(jié)合腔中的S1，S2，S4和S1′結(jié)合位點(diǎn)緊密結(jié)合時(shí)，化合物將有很高的酶抑制活性.為了克服肽類抑制劑的不良藥代性質(zhì)，Kenichi Akaji課題組通過分子模擬，以絲氨酸及苯基異絲氨酸為骨架（圖1），設(shè)計(jì)了一系列具有良好活性的非肽類小分子抑制劑[14-15]，值得進(jìn)一步深入研究.

圖1 絲氨酸及苯基異絲氨酸類抑制劑與3CL蛋白酶結(jié)合模式Fig.1 Binding mode of serines and phenylisoserines inhibitors with 3CL protease

定量構(gòu)效關(guān)系（QSAR）研究在現(xiàn)代藥物設(shè)計(jì)中起著重要的作用.二維定量構(gòu)效關(guān)系（2D-QSAR）和三維定量構(gòu)效關(guān)系（3D-QSAR）已經(jīng)得到了廣泛應(yīng)用[16-19]，但這兩種方法均有其自身局限性.傳統(tǒng)的2D-QSAR在分析過程中涉及大量的物化參數(shù)，需對(duì)這些參數(shù)進(jìn)行選擇及組合，過程復(fù)雜[20].3DQSAR中如比較分子力場(chǎng)分析（CoMFA）方法，雖然包含了更為明確的物理化學(xué)信息，同時(shí)也能夠間接反映分子與靶點(diǎn)之間的非鍵相互作用特征[21]，但其模建過程依賴于分子的疊合，大大增加了工作量.1997年，Hurst和Heritage建立了分子全息定量構(gòu)效關(guān)系(HQSAR)方法[22]，無需疊合分子構(gòu)象，計(jì)算簡(jiǎn)單、快捷，所建模型具有較高的預(yù)測(cè)能力，已被廣泛用于構(gòu)效關(guān)系研究中[23-26].

本文采用HQSAR方法研究了絲氨酸及苯基異絲氨酸類3CL蛋白酶小分子抑制劑的構(gòu)效關(guān)系，建立了具有較好預(yù)測(cè)能力的構(gòu)效關(guān)系模型.并以此為指導(dǎo)，對(duì)該類抑制劑進(jìn)一步結(jié)構(gòu)優(yōu)化，經(jīng)模型預(yù)測(cè)，所設(shè)計(jì)化合物對(duì)3CL蛋白酶均有較高的抑制活性. 本研究為此類3CL蛋白酶小分子抑制劑的設(shè)計(jì)與優(yōu)化提供了指導(dǎo)，也為新型SARS-CoV-2 3CL蛋白酶抑制劑的設(shè)計(jì)提供參考.

1 計(jì)算方法

1.1 化合物及活性數(shù)據(jù)用于HQSAR分析的絲氨酸及苯基異絲氨酸類SARS-CoV 3CL蛋白酶抑制劑源自Kenichi Akaji課題組[14-15].分子公共結(jié)構(gòu)見圖2，其中化合物37，38，39，40在表1中額外列出.所選取的40個(gè)化合物的結(jié)構(gòu)見表1，活性數(shù)據(jù)見表2.

表1 40種化合物結(jié)構(gòu)Tab.1 Structure of 40 compounds

表2 39個(gè)化合物生物活性的實(shí)驗(yàn)值、計(jì)算值和誤差Tab.2 Biocative experiment data，calculated value and deviatoin of 39 compounds

圖2 絲氨酸及苯基異絲氨酸類3CL蛋白酶抑制劑公共結(jié)構(gòu)（1～36）Fig.2 The structure of serine and phenylisoserine 3CL protease inhibitors(1～36)

活性指標(biāo)用抑制3CL蛋白酶活性的半數(shù)有效濃度IC50的負(fù)對(duì)數(shù)pIC50(-lgIC50)表示.在初步構(gòu)效關(guān)系研究中，化合物2的R為氫，異于其它化合物的烴基取代基，在HQSAR模型中為異常點(diǎn)，予以去除.故本文最終選取了30個(gè)化合物作為訓(xùn)練集構(gòu)建模型，隨機(jī)選取9個(gè)化合物作為測(cè)試集（表1中帶*的分子），用于檢驗(yàn)?zāi)Ｐ皖A(yù)測(cè)能力.

1.2 分子建模本文所有計(jì)算工作均在Tripos公司的SYBYL7.2分子模擬軟件上進(jìn)行.計(jì)算中各項(xiàng)參數(shù)除特別指明外均采用缺省值.首先使用SYBYL/Base Builder模塊構(gòu)建化合物結(jié)構(gòu)，隨后利用Minimize模塊中的Tripos分子力場(chǎng)進(jìn)行能量?jī)?yōu)化，采用Gasteiger-Huckel方法計(jì)算原子電荷，能量收斂標(biāo)準(zhǔn)為0.0209 kJ/mol，迭代次數(shù)為5 000次，優(yōu)化后得到的各分子最低能量構(gòu)象用于HQSAR分析.

1.3 HQSAR模型的構(gòu)建采用SYBYL/HQSAR模塊產(chǎn)生分子全息.首先，將訓(xùn)練集中的各分子剪切為包含一定原子個(gè)數(shù)的分子碎片，碎片的結(jié)構(gòu)類型即為碎片區(qū)分參數(shù)，SYBYL中提供的碎片區(qū)分參數(shù)為原子類型(A)、化學(xué)鍵類型(B)、連接性(Co)、氫原子(H)、手性(Ch)及氫鍵供體或受體(DA).碎片中的原子數(shù)目即為碎片大小參數(shù)(fragment size parameters)，由碎片中最小原子數(shù)N和最大原子數(shù)M設(shè)定.一般而言，碎片長(zhǎng)度在1～3為較小分子碎片，4～7為中等原子碎片，8～10為較大原子碎片[27].隨后，基于循環(huán)冗余校驗(yàn)算法將產(chǎn)生的分子碎片映射為0～231的偽隨機(jī)整數(shù)，使得每個(gè)分子可用一定長(zhǎng)度的整數(shù)串進(jìn)行表示.繼而基于哈希法將其轉(zhuǎn)換為具有相同長(zhǎng)度的整數(shù)串，即分子全息[28]. 在HQSAR模塊中，提供了12個(gè)缺省質(zhì)數(shù)(53，59，61，71，83，97，151，199，257，307，353和401)作為全息長(zhǎng)度[29].通過調(diào)整不同的碎片大小、全息長(zhǎng)度等參數(shù)可以得到不同的分子全息.隨后，將各個(gè)分子的分子全息作為自變量，活性作為因變量，采用偏最小二乘法(PLS)進(jìn)行線性回歸分析.采用留一法(LOO)進(jìn)行交叉驗(yàn)證，通過交叉驗(yàn)證相關(guān)系數(shù)q2確定最佳主成分?jǐn)?shù)，選擇q2最大時(shí)的主成分?jǐn)?shù)為最佳主成分?jǐn)?shù)(NOC)，應(yīng)用最佳主成分?jǐn)?shù)進(jìn)行非交叉驗(yàn)證PLS分析得到HQSAR模型.

2 結(jié)果與討論

2.1 碎片區(qū)分參數(shù)的確定由于不同的碎片區(qū)分參數(shù)及其組合所包含的信息不同，因此確定適宜的碎片區(qū)分參數(shù)對(duì)HQSAR模型至關(guān)重要.首先，在默認(rèn)長(zhǎng)度(4～7)下及較大的分子碎片(1～9)下考察不同碎片區(qū)分參數(shù)以及分子全息長(zhǎng)度對(duì)模型的影響，結(jié)果見表3.由于本文所研究的化合物具有手性，因此碎片區(qū)分參數(shù)中的手性(Ch)為必選項(xiàng).由表3的數(shù)據(jù)可知，當(dāng)以原子類型(A)、化學(xué)鍵類型(B)、連接性(Co)和手性(Ch)作為碎片區(qū)分參數(shù)時(shí)，可以得到較好的模型.此時(shí)，留一法交叉驗(yàn)證相關(guān)系數(shù)q2為0.504，非交叉驗(yàn)證相關(guān)系數(shù)r2=0.937，最佳主成分?jǐn)?shù)(NOC)為6，最佳分子全息長(zhǎng)度(BL)為97.

表3 碎片區(qū)分參數(shù)及全息長(zhǎng)度對(duì)模型的影響Tab.3 The influence of various fragment distinction parameters and hologram length on HQSAR

2.2 碎片大小的確定由于不同大小的分子碎片所包含的分子結(jié)構(gòu)信息有所區(qū)別，所以分子碎片大小也會(huì)對(duì)HQSAR模型造成影響.為得到最優(yōu)模型，以A/B/Co/Ch為最優(yōu)碎片區(qū)分參數(shù)，以不同大小的分子碎片建立HQSAR模型，表4表明了不同碎片長(zhǎng)度對(duì)HQSAR結(jié)果的影響.當(dāng)分子碎片為1～15時(shí)可得最優(yōu)模型，q2達(dá)到0.604，r2為0.904.

表4 不同碎片大小對(duì)HQSAR分析結(jié)果的影響Tab.4 Influences of various fragment parameters for HQSAR

2.3 最佳模型的評(píng)價(jià)與檢驗(yàn)在QSAR分析中，r2和q2是衡量模型效果的主要標(biāo)準(zhǔn).通常當(dāng)r2≥0.9時(shí)認(rèn)為該模型具有良好的擬合能力，當(dāng)q2≥0.5時(shí)則表明模型具有良好的預(yù)測(cè)能力.本文所得最佳模型的r2為0.904，q2為0.604，表明所建模型具有良好的相關(guān)性和較強(qiáng)的預(yù)測(cè)能力.此外，模型的非交叉驗(yàn)證標(biāo)準(zhǔn)誤差(SEE)也是衡量模型的重要標(biāo)準(zhǔn)，本模型SEE為0.125，活性數(shù)據(jù)的數(shù)值范圍約為1.57，比值約為7.96%，小于10%的標(biāo)準(zhǔn)，故所建模型良好[30].為進(jìn)一步驗(yàn)證模型的預(yù)測(cè)能力及穩(wěn)定性，采用最優(yōu)模型對(duì)測(cè)試集進(jìn)行了活性預(yù)測(cè)，結(jié)果見表2和圖3.由計(jì)算結(jié)果可知，預(yù)測(cè)值與實(shí)驗(yàn)值偏差較小，其預(yù)測(cè)相關(guān)系數(shù)r2pred為0.723，表明此模型具有良好的預(yù)測(cè)能力和穩(wěn)定性，可用于此類化合物抗SARS-CoV活性預(yù)測(cè).

圖3 訓(xùn)練集和測(cè)試集樣本的實(shí)際值與預(yù)測(cè)值相關(guān)關(guān)系Fig.3 Comparison of the experiment value and the predicted value of the training set and test set

綜上所述，當(dāng)采用碎片區(qū)分參數(shù)為A/B/Co/Ch，碎片大小為1～15，最佳全息長(zhǎng)度為71，最佳組分?jǐn)?shù)為4時(shí)可得最佳HQSAR模型，此時(shí)預(yù)測(cè)相關(guān)系數(shù)為0.723，交叉驗(yàn)證相關(guān)系數(shù)q2為0.604，非交叉驗(yàn)證相關(guān)系數(shù)r2為0.904，標(biāo)準(zhǔn)誤差SEE為0.125.

2.4 HQSAR模型圖解釋由HQSAR模型得到的色碼圖，采用不同顏色表征分子中不同原子或片段對(duì)活性的貢獻(xiàn)，從而能直觀地指導(dǎo)分子結(jié)構(gòu)的改造與優(yōu)化.部分化合物HQSAR活性色碼圖如圖4中色譜帶所示，從紅色到綠色的變化，表示基團(tuán)或原子對(duì)分子活性的貢獻(xiàn)由不利到有利. 紅色、橙紅色表示該部位為分子活性提供負(fù)貢獻(xiàn)；白色代表該部位對(duì)分子活性貢獻(xiàn)中等；黃色、綠色表示該部位原子或基團(tuán)對(duì)整個(gè)分子活性提供正貢獻(xiàn).

圖4 HQSAR活性貢獻(xiàn)圖Fig.4 HQSAR contribution maps

在a類化合物（1～15）中，基團(tuán)R作用于酶活性中心的S1′結(jié)合位點(diǎn).當(dāng)R為苯甲?；鶗r(shí)得到化合物1(pIC50=3.903)，由其色碼圖可以看出，其分子主鏈呈黃色及貢獻(xiàn)度中等的白色.而當(dāng)R為苯乙酰基，叔丁?；?-甲?；拎r(shí)，得到的化合物4，7，10的pIC50分別為3.187，3.347和3.745，活性大大降低，其色碼圖中顯示出不利于活性的橙紅色，表明以苯甲酰基作為R取代基，有利于提高化合物的活性.

對(duì)化合物1的苯甲?；M(jìn)行修飾，當(dāng)R基團(tuán)為2-甲基-6-硝基苯甲?；?，得到化合物15(pIC50=4.187)，其活性明顯提高，且分子主鏈存在呈正貢獻(xiàn)的黃色及綠色結(jié)構(gòu).而當(dāng)R為2，4，6-三氯苯甲?；?，5-二硝基苯甲?；鶗r(shí)，得到化合物6和9，其pIC50分別為3.260和3.602，和化合物1相比，其活性均出現(xiàn)了不同程度的下降，且色碼圖中出現(xiàn)了呈負(fù)貢獻(xiàn)的區(qū)域.因此，以2-甲基-6-硝基苯甲?；鳛镽取代基，有利于化合物活性的提高.

在b類化合物（16～24）中，基團(tuán)R則作用于酶活性中心的S4結(jié)合位點(diǎn).當(dāng)R為單甲氧基或雙甲氧基取代肉桂酰基時(shí)，得到化合物17，19，20，21和22，其pIC50分別為3.620，3.810，3.812，4.000，4.009.由色碼圖可以看出，在R基末端苯環(huán)3位引入甲氧基后，其色碼圖中呈正貢獻(xiàn)的黃綠色部分均比引入前有所增加，活性也有所提高；當(dāng)R為3，4，5-三甲氧基肉桂?；鶗r(shí)得到化合物24(pIC50=4.131)，其活性高于其它b類分子，且色碼圖中呈正貢獻(xiàn)的黃綠色區(qū)域也遠(yuǎn)多于其它b系化合物，表明3，4，5-三甲氧基肉桂?；芡?CL蛋白酶活性中心的S4結(jié)合位點(diǎn)形成有利于活性的相互作用.

與具有絲氨酸骨架的a，b兩系列分子不同的是，c類分子（25～30）為苯基異絲氨酸類化合物，不僅能同3CL蛋白酶活性中心的S1′和S4結(jié)合位點(diǎn)發(fā)生作用，其苯基還能同S2結(jié)合位點(diǎn)形成相互作用.比較c類化合物25和a類化合物1的色碼圖可以明顯看出，當(dāng)與酶活性中心S1′和S4結(jié)合位點(diǎn)作用的取代基相同的情況下，化合物25的色碼圖中對(duì)活性有利的黃色及綠色部分遠(yuǎn)遠(yuǎn)多于化合物1的正貢獻(xiàn)部分，同時(shí)，化合物25（pIC50=4.367）的活性較化合物1（pIC50=3.903）提高了一個(gè)數(shù)量級(jí)，說明以苯基異絲氨酸為骨架進(jìn)行分子結(jié)構(gòu)改造與優(yōu)化，能更有效地提高化合物的酶抑制活性.此外，化合物25是c類化合物中活性最好的，其色碼圖中對(duì)活性有利的黃綠色部分也多于其它c(diǎn)類分子，表明在c類化合物中苯甲?；廊荒芡富钚灾行牡腟1′位點(diǎn)形成有利于活性的相互作用.

綜上所述，以苯基異絲氨酸為骨架，苯甲?；蛉〈郊柞；鳛楹?CL蛋白酶活性中心S1′位點(diǎn)作用的取代基，而與S4位點(diǎn)作用的取代基優(yōu)選取代肉桂?；?，均能有效提高化合物的活性，這為該類化合物的進(jìn)一步結(jié)構(gòu)改造與優(yōu)化指明了方向.

2.5 基于HQSAR結(jié)論的分子設(shè)計(jì)基于以上色碼圖分析，我們以苯基異絲氨酸為骨架進(jìn)行結(jié)構(gòu)修飾，其中，和S1′結(jié)合位點(diǎn)作用的R1基團(tuán)分別為苯甲酰基，2-甲基-6-硝基苯甲酰基及肉桂酰基，和S4結(jié)合位點(diǎn)作用的R2基團(tuán)則為肉桂?；凹籽趸〈夤瘐；?，并采用所建立的最佳HQSAR模型預(yù)測(cè)所設(shè)計(jì)分子的活性，分子結(jié)構(gòu)及預(yù)測(cè)值見表5和圖5.由可以看出，基于所建模型設(shè)計(jì)的化合物均具有較優(yōu)的預(yù)測(cè)活性，且5個(gè)化合物預(yù)測(cè)活性均大于等于現(xiàn)有化合物測(cè)試活性.以上結(jié)果支持了以苯基異絲氨酸為骨架進(jìn)行分子結(jié)構(gòu)改造，能更有效地提高化合物酶抑制活性的理論假設(shè).

表5 基于最優(yōu)HQSAR結(jié)果所設(shè)計(jì)分子的結(jié)構(gòu)及活性預(yù)測(cè)值Tab.5 Structure and calculated values of the designed molecules

圖5 基于HQSAR設(shè)計(jì)的化合物通用結(jié)構(gòu)Fig.5 General structure of compounds designed by HQSAR

3 結(jié)論

本文采用HQSAR方法對(duì)絲氨酸及苯基異絲氨酸類SARS-CoV 3CL蛋白酶抑制劑進(jìn)行了分子全息定量構(gòu)效關(guān)系研究，得到的最優(yōu)HQSAR模型交叉驗(yàn)證相關(guān)系數(shù)q2為0.604，非交叉驗(yàn)證相關(guān)系數(shù)r2為0.904，標(biāo)準(zhǔn)偏差SEE為0.125.對(duì)由9個(gè)化合物構(gòu)成的測(cè)試集進(jìn)行了預(yù)測(cè)，測(cè)試集的預(yù)測(cè)相關(guān)系數(shù)r2pred為0.723，表明該模型具有較高的預(yù)測(cè)能力和擬合能力.通過HQSAR色碼圖直觀展示的不同取代基對(duì)化合物活性呈現(xiàn)的有利及不利貢獻(xiàn)，設(shè)計(jì)出一組具有較高抗3CL蛋白酶活性的苯基異絲氨酸類抑制劑，為此類化合物的進(jìn)一步結(jié)構(gòu)改造與優(yōu)化提供指導(dǎo)，也為以3CL蛋白酶為潛在靶標(biāo)的SARS-CoV-2冠狀病毒抑制劑的設(shè)計(jì)與結(jié)構(gòu)優(yōu)化提供參考.