基于計算機輔助水解的中藥大豆寡肽的ETA拮抗活性預測

喬連生+蔣蘆荻+雒剛剛+路芳+陳艷昆+王靈芝+李貢宇+張燕玲

[摘要] 中藥寡肽是中藥發(fā)揮藥效的關鍵組分之一，系統(tǒng)地研究中藥寡肽的組成及其藥效是中藥物質基礎及作用機制研究的關鍵。該研究擬基于計算機輔助水解和分子對接等模擬技術，以中藥大豆為研究載體，解析其降壓寡肽成分，并預測其潛在的內皮素受體A（endothelin receptor A，ETA）拮抗活性。該文通過收集大豆中的儲存蛋白序列，基于計算機輔助水解方法進行蛋白的模擬消化水解，并將水解獲得的寡肽構建虛擬結構數(shù)據(jù)庫。隨后該研究構建了ETA肽類拮抗劑的藥效團模型，包括1個疏水特征、1個負電中心、1個芳環(huán)特征和5個排除體積。同時，利用同源模建方法構建ETA的蛋白三維結構模型，并將其用于分子對接研究，采用一致性打分評價化合物的ETA拮抗活性。通過構建ETA肽類拮抗劑的藥效團模型、同源模建的三維蛋白結構以及分子對接模型，共預測獲得27條可能具有潛在ETA拮抗活性的寡肽分子。該文進一步分析關鍵氨基酸GLN165，并結合文獻，說明其對拮抗劑活性的重要性。計算機輔助水解方法可以高效地對已知序列結構的中藥蛋白進行模擬水解，結合分子模擬模型，能進一步辨識中藥寡肽潛在的生物學活性。該研究為快速、高效開展中藥來源的肽類物質的活性機制研究提供了方法學依據(jù)。

[關鍵詞] 計算機輔助水解； 寡肽； 大豆； 內皮素受體A； 藥效團； 同源模建； 分子對接

[Abstract] Oligopeptides are one of the the key pharmaceutical effective constituents of traditional Chinese medicine（TCM）. Systematic study on composition and efficacy of TCM oligopeptides is essential for the analysis of material basis and mechanism of TCM. In this study， the potential anti-hypertensive oligopeptides from Glycine max and their endothelin receptor A （ETA） antagonistic activity were discovered and predicted based on in silico technologies.Main protein sequences of G. max were collected and oligopeptides were obtained using in silico gastrointestinal tract proteolysis. Then， the pharmacophore of ETA antagonistic peptides was constructed and included one hydrophobic feature， one ionizable negative feature， one ring aromatic feature and five excluded volumes. Meanwhile， three-dimensional structure of ETA was developed by homology modeling methods for further docking studies. According to docking analysis and consensus score， the key amino acid of GLN165 was identified for ETA antagonistic activity. And 27 oligopeptides from G. max were predicted as the potential ETA antagonists by pharmacophore and docking studies.In silico proteolysis could be used to analyze the protein sequences from TCM. According to combination of in silico proteolysis and molecular simulation， the biological activities of oligopeptides could be predicted rapidly based on the known TCM protein sequence. It might provide the methodology basis for rapidly and efficiently implementing the mechanism analysis of TCM oligopeptides.

[Key words] in silico proteolysis； oligopeptides； Glycine max； ETA； pharmacophore； homology modeling； docking

中藥蛋白是中藥中一類常見的活性物質，具有明顯而特殊的藥理活性[1]。從兩千多年前阿膠的應用，到水蛭素的制備和大豆多肽的上市，中藥蛋白和生物活性肽的應用具有悠久的歷史，已有大量的藥性、藥理和臨床數(shù)據(jù)。而中藥蛋白具有物理化學性質不穩(wěn)定、易水解和易氧化的特點，因此，基于傳統(tǒng)的中藥化學方法，通過對中藥蛋白提取分離酶解獲取中藥活性肽，往往具有分離困難、耗時費力的特點。而計算機輔助水解是基于生物信息學方法，預測蛋白質序列在特定蛋白酶或化學物質在給定條件下的切割位點的技術，該技術具有快速高效可靠的寡肽預測能力[2]，已較為廣泛的應用于食品工業(yè)研究[3-4]。

中藥大豆是豆科植物大豆Glycine max的成熟種子，現(xiàn)代研究表明大豆活性肽具有良好的降壓活性[5]。內皮素受體A（endothelin receptor A，ETA）是治療高血壓和肺動脈高壓的經典靶標，其第一代拮抗劑即為肽類化合物[6]。因此，基于大豆寡肽開發(fā)ETA拮抗劑具有較高的研究價值。本研究擬利用計算機輔助水解技術，進行中藥大豆的模擬消化水解，并結合基于藥效團和分子對接的虛擬篩選技術，預測中藥大豆寡肽的ETA拮抗活性。

1 材料與方法

1.1 大豆蛋白的模擬水解 本研究在Uniport（http：//www.uniprot.org/）中檢索已解析的大豆種子儲存蛋白的多肽序列，共獲得5條大豆球蛋白（glycinin）和3條伴大豆球蛋白（conglycinin）序列，主要包括大豆球蛋白G，大豆球蛋白G1-G4，β-伴大豆球蛋白的α，α′和β鏈。其為大豆中的主要儲存蛋白，占分離蛋白的80%以上[7]。基于獲得的蛋白序列，本研究采用BIOPEP（http：//www.uwm.edu.pl/biochemia/index.php/pl/biopep）進行計算機輔助模擬水解，以預測蛋白序列的切割位點。本研究采用模擬胃腸道的水解策略，選擇胃蛋白酶（pH 1.3）、胰蛋白酶和糜蛋白酶作為消化酶[8]。針對獲得的寡肽序列，利用Discovery Studio 4.0（DS）進行寡肽三維結構數(shù)據(jù)庫的構建。

1.2 ETA肽類拮抗劑的藥效團模型的構建 為了更可靠的篩選ETA的肽類拮抗劑，本研究擬采用DS中的3D-QSAR pharmacophore（Hypogen）模塊，構建ETA肽類拮抗劑的藥效團模型。構建藥效團的訓練集和測試集中的50個活性化合物均來源于文獻報道[9]，其包括直鏈肽和環(huán)肽2類化學結構，其活性值分布5個數(shù)量級。

對訓練集中20個化合物在CHARMm中進行三維構象生成和能量最小化處理。隨后，對其進行藥效特征分析，發(fā)現(xiàn)活性化合物中包括的藥效特征有氫鍵受體（hydrogen bond acceptor，A）、親脂性氫鍵受體（hydrogen bond acceptor lipid，Ali）、氫鍵供體（hydrogen bond donor，D）、疏水基團（hydrophobic region，H）、脂性疏水基團（hydrophobic aliphatic region，Hal）、芳香疏水基團（hydrophobic aromatic region，Har）、芳香環(huán)（ring aromatic，R）、負電基團（ionizable negative group，N）。

基于特征分析，進一步構建藥效團模型。首先，對訓練集中活性化合物進行構象分析，采用Best模式，設置能量收斂值為20 kcal·mol-1，最大構象數(shù)目255個。選擇生成的5個藥效特征為氫鍵受體、氫鍵供體、疏水基團、芳香環(huán)和負電基團，同時，模型中通過增加5個排除體積（maximum excluded volumes，Ev）以提高模型的特異性，并設置化合物的活性不確定度為1.5，利用Hypogen模塊構建10個候選藥效團模型。

針對生成的藥效團模型需進行進一步的模型評價。利用剩余的30個活性化合物及從Binding Database（http：//www.bindingdb.org/bind/index.jsp）中隨機挑選的90個非活性化合物，以1∶3的比例生成測試集，并對測試集分子進行同訓練集分子相同的計算處理。隨后，利用生成的候選藥效團模型與測試集化合物進行匹配篩選，計算相應的有效命中率（hit rate of active compounds，HRA）、辨識有效性指數(shù)（identify effective index，IEI）、綜合評價指數(shù)（comprehensive appraisal index，CAI）等指標對藥效團優(yōu)劣進行評價[10]。3個外部評價指標計算公式如下所示（1）～（3），其中TD代表測試集中化合物總數(shù)，TA代表測試集中活性化合物數(shù)目。Ha代表命中的活性化合物數(shù)目，Ht代表命中化合物總數(shù)。本文綜合利用藥效團模型的多種評價參數(shù)，選出最優(yōu)的藥效團模型。利用最優(yōu)的藥效團模型進行大豆寡肽庫的虛擬篩選，其中搜索模式為Best Search、匹配模式為柔性匹配，初步獲得潛在的活性寡肽。

HRA=Ha/TA×100%（1）

IEI=（Ha/Ht）/（TA/TD）（2）

CAI=IEI×HRA（3）

1.3 ETA的同源模建及分子對接 ETA（Uniprot/P25101）由427個氨基酸殘基組成，是典型的G蛋白偶聯(lián)受體（G protein-coupled receptor，GPCR）中A家族的7次跨膜蛋白。通過對已解析的同家族蛋白晶體進行序列比對，分析其同源性和相似性的高低，選擇相似性較高的蛋白作為模板蛋白。本研究同源模建使用Modeller 9.13程序，采用多模板建模方法，生成5個候選模型進行進一步優(yōu)化。首先，基于逐級能量最小化方法，先約束所有重原子，優(yōu)化氫原子，再約束保守跨膜區(qū)，優(yōu)化側鏈及Loop區(qū)，最后優(yōu)化跨膜區(qū)。能量最小化過程采用最陡下降法和共軛梯度法進行優(yōu)化。最后基于拉氏圖和ERRAT分析，利用SAVES服務器（http：//services.mbi.ucla.edu/SAVES/）對同源模建蛋白進行評價。

結合對蛋白空腔的自動識別，同時參考GPCR同家族活性口袋的位置，定義ETA蛋白的活性口袋。由于ETA活性口袋中可能存在多種結合模式[11]，為保證分子對接結果的準確性，本研究選擇較大的活性口袋半徑。本研究利用DS中的Libdock模塊進行分子對接。首先，以陽性藥作為測試，分別對接進入活性口袋中，并同時計算5種打分函數(shù)，采用一致性打分的方法，評價對接模型的可靠性，并分析ETA肽類拮抗劑可能存在的作用模式。隨后，利用上述的對接模型進一步篩選藥效團的命中結果，篩選潛在的ETA肽類拮抗劑。

2 結果與討論

2.1 大豆寡肽庫的構建 通過模擬水解，本研究共獲得449條多肽，其中各段寡肽的數(shù)目分布見圖1，由于2～6肽可能具有較好的相對分子質量與成藥性[2]，所以選擇2～6肽共299條作為本研究的主要研究對象。通過對模擬水解結果與已知的活性肽數(shù)據(jù)庫BIOPEP比對，共發(fā)現(xiàn)已知活性肽32條，其主要的生物學活性包括ACE抑制作用和抗氧化作用等，其結果見表1，表明大豆寡肽具有較為廣泛的藥理活性。

2.2 ETA肽類拮抗劑藥效團模型的建立與評價 基于Hypogen的藥效團建模方法主要包括2類內部評價參數(shù)，即訓練集活性化合物的預測結果與真實活性的相關系數(shù)Corrtrain，以及費用函數(shù)（Cost）對所構建的藥效團模型的顯著性評價。本研究構建的10個藥效團內部評價參數(shù)見表2。模型的Null Cost為425.75，ΔCost為Null Cost與Total Cost的差值，ΔCost大于60說明模型有90%以上的顯著性。

利用測試集對10個候選藥效團進行評價，計算相應外部評價參數(shù)，包括藥效團的活性化合物命中率HRA、辨識有效性指數(shù)IEI和綜合評價指數(shù)CAI以及測試集中活性化合物的預測值與真實值之間的相關系數(shù)Corrtest，其結果見表2。本研究生成的藥效團模型主要包括2類，即HNR型和DHR型。HNR型主要包括1個疏水特征，1個負電中心和1個芳環(huán)基團，而DHR型主要包括1個氫鍵受體，1個疏水特征和1個芳香基團。通過對比其CAI值差異，發(fā)現(xiàn)HNR型具有更好的測試集篩選能力。進一步分析藥效團對測試集的活性預測能力，發(fā)現(xiàn)1號和5號藥效團具有相當?shù)幕钚曰衔锉孀R能力。通過比較2個藥效團的空間結構差異，發(fā)現(xiàn)2個藥效團中HNR 3個特征之間的空間距離及角度具有相似性，同時，5號藥效團的D特征與其排除體積相重疊，因此最終選擇1號藥效團作為最優(yōu)藥效團對大豆寡肽庫進行虛擬篩選，結果共獲得64條寡肽。其中，VVF與藥效團具有較好的匹配，其匹配值為0.85，見圖2。

2.3 ETA同源模建及分子對接研究 通過同源比對，選擇3個蛋白晶體結構作為模板進行同源模建，包括Delta Opioid Receptor（PDB：4N6H_A），Neurotensin Receptor 1（PDB：4BV0_A）和orphanin FQ receptor（PDB：4EA3_A）。3個受體均來源于GPCR的A家族，且其與ETA的序列相似度 （similarity）分別為40.0%，39.8%，38.8%，其一致性（identity）分別為26%，24%，25%。利用獲得的3個模板進行多模板建模，通過能量最小化處理，利用拉氏圖和ERRAT評價模型的優(yōu)劣，最終模型和評價參數(shù)見圖3。其中，拉氏圖的最適殘基量為96.96%，ERRAT的值為91.41，顯示模建結構具有良好的立體化學參數(shù)和空間結構。

ETA分子對接結果：本研究通過分析ETA的空腔結構，自動識別了8個活性口袋，通過文獻中關鍵氨基酸的描述和參考GPCR同家族的位點情況，最終確定對接口袋，其半徑為15 。先將訓練集中20個陽性寡肽與活性口袋對接，20個化合物均能對接進入活性口袋中，進一步計算Jain，-PLP2，-PMF，-PMF04和-PLP1 5種打分函數(shù)，綜合已有的LibDockScore打分，進行一致性打分評價[12]。其中，一致性打分大于3的化合物為95.0%，說明了對接模型和一致性打分的可靠性[13]。隨后，分析ETA與陽性藥發(fā)生相互作用的高頻氨基酸，產生氫鍵作用的氨基酸主要包括：LYS329，ARG326，ASP351和GLN165，產生疏水作用的氨基酸主要有LEU347，LYS140，LEU141，ARG145，ARG326，LYS329，與文獻報道基本一致[14-15]。

進一步將測試集中30個活性分子與活性口袋對接，一致性打分大于3的化合物為96.7%，說明對接模型和一致性打分評價對化合物具有較好的預測能力。發(fā)生氫鍵作用的關鍵氨基酸主要有GLN165，LYS329和ASP351，產生疏水作用的高頻氨基酸有ARG326，LEU322和LYS329等，與文獻報道及訓練集化合物基本一致。其中高活性化合物Com1385與ETA之間的相互作用見圖4。

利用構建的分子對接模型和一致性打分方法，篩選被藥效團命中的64條寡肽。27條寡肽的一致性打分為6，提示大豆虛擬水解獲得的寡肽可能具有較好的ETA拮抗活性。其中，VVF能與關鍵氨基酸具有良好的匹配，VVF中苯丙氨酸的羧基與藥效團的負電中心匹配，而其羰基氧與GLN165發(fā)生氫鍵相互作用，而其羥基氫能與ASP351發(fā)生氫鍵相互作用，因此，VVF的苯丙氨酸羧基可能是發(fā)揮活性的重要基團[16]，其結果見圖5。

3 結論

本研究基于計算機輔助水解方法，對大豆中的主要蛋白進行模擬消化水解，獲得2～6肽的寡肽基本結構共299條。在此基礎上，通過構建ETA肽類拮抗劑的藥效團，以及ETA蛋白同源模建與分子對接，篩選并預測大豆寡肽的ETA拮抗活性，共獲得27條可能具有ETA拮抗活性的寡肽分子，其結果有待進一步的活性驗證。本文通過基于生物信息學和計算機輔助藥物設計的方法，從大豆蛋白中篩選潛在的ETA拮抗活性寡肽，為探討中藥蛋白發(fā)揮藥效的作用機制探討提供了一條新的研究方法，為開展中藥蛋白的酶解分離或藥理驗證提供指導。

[參考文獻]

[1] 鄒吉利，徐南平.中藥活性多肽研究進展[J].湖北中醫(yī)藥大學學報， 2012， 14（4）：66.

[2] 鄒平. 基于生物信息學與QSAR及分子對接的菜粕活性肽篩選及活性研究[D]. 杭州：浙江大學， 2014.

[3] Dziuba B， Dziuba M. New milk protein-derived peptides with potential antimicrobial activity： an approach based on bioinformatic studies[J]. Int J Mol Scis， 2014， 15（8）：14531.

[4] 孫國威，樂國偉，施用暉.模擬酶解大豆7S、11S蛋白及其抗氧化活性的研究[J].食品工業(yè)科技， 2010（7）：101.

[5] 于婷婷，韓飛，陳光.大豆降壓肽研究進展[J].糧油食品科技， 2008， 16（2）：27.

[6] 楊菁.內皮素受體拮抗劑的設計與合成[D].北京：中國人民解放軍軍事醫(yī)學科學院，2008.

[7] 田琨，管娟，邵正中，等.大豆分離蛋白結構與性能[J].化學進展， 2008， 20（4）：565.

[8] Lafarga T， O′Connor P， Hayes M. Identification of novel dipeptidyl peptidase-IV and angiotensin-I-converting enzyme inhibitory peptides from meat proteins using in silico analysis[J]. Peptides， 2014， 59（9）：53.

[9] Fukami T， Yamakawa T， Niiyama K， et al. Synthesis and structure-activity relationships of 2-substituted D-tryptophan-containing peptidic endothelin receptor antagonists： importance of the C-2 substituent of the D-tryptophan residue for endothelin A and B receptor subtype selectivity[J]. J Med Chem， 1996， 39（12）：2313.

[10] Jiang L， Zhang X， Chen X， et al. Virtual screening and molecular dynamics study of potential negative allosteric modulators of mGluR1 from Chinese herbs[J]. Molecules， 2015， 20：12769.

[11] Mey J G R D， Compeer M G， Pieter L， et al. ETA -receptor antagonists or allosteric modulators[J]. Trends Pharmacol Sci， 2011， 32（6）：345.

[12] 劉維國.抗藥物依賴功效藥效團模型的構建及其虛擬篩選[D].北京：中央民族大學， 2013.

[13] Krovat E M， Langer T. Impact of scoring functions on enrichment in docking-based virtual screening： an application study on renin inhibitors[J]. J Chem Inf Comp Sci， 2004， 44（3）： 1123.

[14] Breu V， Hashido K， Broger C， et al. Separable binding sites for the natural agonist endothelin-1 and the non-peptide antagonist bosentan on human endothelin-A receptors[J]. Eur J Biochem， 1995， 231（1）：266.

[15] Jens L， Alexander O， Michael B， et al. Structural determinants for selective recognition of peptide ligands for endothelin receptor subtypes ETA and ETB[J]. J Pept Sci， 2009， 15（7）：479.

[16] Khuraijam Dhanachandra S， Karthikeyan M. Molecular modeling， quantum polarized ligand docking and structure-based 3D-QSAR analysis of the imidazole series as dual AT1 and ETA receptor antagonists[J]. Acta Pharm Sin， 2013， 34（12）：1592.