基于埃博拉病毒糖蛋白單體結(jié)構(gòu)從中藥天然產(chǎn)物庫中篩選新型抗病毒化合物

史海龍，黃 月，程 怡，王月雯，馮雪松，晁 旭

陜西中醫(yī)藥大學，西咸新區(qū) 712046

埃博拉病毒(ebola virus，EBOV)是一種包膜絲狀病毒,其引起的急性出血熱是多發(fā)于靈長類動物的人畜共患傳染病，其病死率高達50%～90%。近幾年非洲西部多次暴發(fā)嚴重的埃博拉疫情，但是迄今為止尚無特效的抗埃博拉出血熱治療藥物獲準上市。EBOV進入宿主細胞是感染的關鍵步驟，由病毒表面跨膜單鏈糖蛋白(single glycoprotein of ebolavirus，EBOV-sGP)介導，這種I 類融合蛋白包含二硫鍵連接的兩個亞基，即表面暴露的GP1與嵌入病毒膜GP2。在通過細胞表面附著因子與宿主細胞結(jié)合后，病毒被內(nèi)吞，內(nèi)體相應的半胱氨酸蛋白酶、組織蛋白酶切割掉GP1亞基的大部分肽鏈，裂解暴露受體結(jié)合位點[1]，這種GP構(gòu)象的轉(zhuǎn)變，有利于宿主-病毒膜融合和隨后病毒基因組的釋放。這種融合機制類似HIV-1病毒和流感病毒[2]。因此，開發(fā)有效干擾EBOV-sGP構(gòu)象轉(zhuǎn)變的抑制劑分子，有助于阻斷病毒入侵和復制。

2019年英國牛津大學結(jié)構(gòu)生物學系Elizabeth E.Fry及其團隊通過使用蛋白表達純化結(jié)晶技術已成功制備EBOV-sGP蛋白，并借助X射線衍射實驗測定其三維結(jié)構(gòu)[3](RCSB Protein Data Bank ID：5JQ7)，提供了計算機虛擬篩選模型。除此之外，由于天然產(chǎn)物，尤其是傳統(tǒng)中藥成分庫，具有結(jié)構(gòu)獨特性與多樣性，已被用于治療各種慢性病和傳染病[4,5]，同樣可用于挖掘抗埃博拉病毒抑制劑。基于以上思考，本研究采用計算機模擬技術，從中藥天然產(chǎn)物庫中篩選EBOV-sGP抑制劑分子。

1 資料與計算方法

分子對接計算均采用Schrodinger 2018的Glide模塊、ADME參數(shù)計算采用Schrodinger 2018的QikProp模塊，MM/GBSA結(jié)合自由能預測[6]采用Schrodinger 2018的Prime MM/GBSA模塊，分子動力學模擬及MM/PBSA結(jié)合自由能預測[7]均采用Gromacs 2019程序包。采用Pymol與Discovery Studio Visualizer軟件繪制受體配體結(jié)合模式圖。

1.1 靶蛋白EBOV-sGP的預處理與活性口袋的確定

從PDB數(shù)據(jù)庫下載EBOV-sGP結(jié)構(gòu)(PDB ID為5JQ7，配體為toremifene)。使用Schrodinger2018軟件的Protein Preparation wizard模塊對5JQ7結(jié)構(gòu)進行預處理。以toremifene的空間坐標為參考，使用Receptor Grid Generation模塊生成Grid文件，用于基于受體結(jié)構(gòu)的虛擬篩選。參考Elizabeth E.Fry課題組發(fā)表的文章[8]，選擇Asn61、Val66、Asp522、Leu554可作為關鍵氨基酸殘基，并且關注Ala101、Leu515、Tyr517、Gly67、Gly102、Arg64、Thr519、Met548、Ile38、Leu186、Leu63與候選藥效分子是否產(chǎn)生疏水相互作用、離子鍵等，并使用Pymol進行可視化處理，如圖1所示。

圖1 EBOV-sGP活性口袋的空間結(jié)構(gòu)

1.2 中藥分子化學結(jié)構(gòu)的預處理

TCMSP數(shù)據(jù)庫[9]中13 144個化合物作為配體庫，采用Schrodinger的LigPrep模塊對各分子進行加氫、能量最小化、幾何優(yōu)化等處理。在力場OPLS3下，計算pH為7.0±2.0范圍的化合物離子狀態(tài)，搜索合理構(gòu)象作為對接起始構(gòu)象，形成可能的互變異構(gòu)體。以同樣方式預處理toremifene。

1.3 基于類藥性評估、分子對接模擬、ADME預測、結(jié)合自由能計算的多輪虛擬篩選

首先，基于類藥性規(guī)則“Lipinski Ro5[10]”，“Verber Ro3[11]”對中藥分子配體庫進行首輪篩選。然后，基于EBOV-sGP活性口袋結(jié)構(gòu)Glide程序?qū)γ蟹肿舆M行第二輪篩選，依次從高通量模式HTVS-docking、標準模式SP-docking、高精度模式XP-docking進行不同精度的分子對接。接著，基于ADME藥代動力學參數(shù)進行第三輪篩選，從分子量(molecular weight，MW)、氫鍵供體(hydrogen bond donor，HB donor)、氫鍵配體(hydrogen bond acceptor，HB acceptor)、溶劑可及表面積(solvent accessible surface area，SASA)、K+通道阻斷相關IC50(IC50value for blockage of HERG K+channels，logHERG)、水溶性(aqueous solubility，logS)、脂水分配系數(shù)(octanol/water partition coefficient，logO/W)、口服生物利用度(oral bioavailability，OB)等評估，篩選具有良好類藥性質(zhì)的化合物。最后，基于Prime MM/GBSA結(jié)合自由能預測進行第四輪篩選，由于高精度分子對接函數(shù)計算仍然無法估算配體-受體復合物結(jié)合的溶劑化效應因素，故使用隱性溶劑模型下MM/GBSA算法優(yōu)化排序。

1.4 分子動力學模擬

為了從動力學和熱力學角度動態(tài)分析EBOV-sGP與命中分子互作模式，使用GROMACS進行分子動力學模擬。蛋白拓撲參數(shù)取自Amber的FF03力場，由PDB2GMX工具獲得；抑制劑的各原子電荷先用Gaussian 09在B3LYP/6-311G**基組水平計算靜電勢，再用Amber中的RESP電荷擬合程序算出。將復合物放入10×10×10 ?3的立方體TIP3P水盒中，并加入Na+和Cl-離子，中和該系統(tǒng)到0.15 mol/L NaCl溶液濃度，隨后進行能量最小化和預平衡，最后在300 K溫度進行50 ns的分子動力學模擬，時間間隔為2 fs，每10 ps保存一次能量和坐標文件。從35～50 ns穩(wěn)定模擬軌跡中提取150幀軌跡文件，計算MM/PBSA結(jié)合自由能。

2 結(jié)果與討論

2.1 虛擬篩選結(jié)果分析

基于靶標EBOV-sGP三維空間結(jié)構(gòu)，從TCMSP數(shù)據(jù)庫13 144個化合物進行多輪虛擬篩選，如圖2所示，其具備有較強的靶標結(jié)合能、良好的成藥性等特質(zhì)。最終確定打分排名前4的化合物MOL006834、MOL000174、MOL002192、MOL012524，具備有較強靶標結(jié)合自由能、良好藥代動力學參數(shù)、較高生物利用度等特質(zhì)，結(jié)果見表1和表2。隨后進行文獻檢索，發(fā)現(xiàn)MOL000174屬于中藥材蒼術的蒼術醇類化合物，是動脈炎病毒的強抑制劑[12]；MOL006834屬于中藥材鴨跖草的阿魏酸酰胺類化合物，鴨跖草具有抑菌、抗氧化、抗病毒、降血脂、鎮(zhèn)痛消炎的作用，其中阿魏酸酰胺類化合物具有抗炎活性[13]；MOL002192屬于中藥材川芎的阿魏酸酯，具有較好的血管擴張作用及抗菌抗癌活性[14]；MOL012524屬于中藥材牛膝的阿魏酰甲氧基酪胺類化合物，在體外藥理實驗中顯示出有效的抗炎活性[15]。從現(xiàn)有的文獻中，4種命中分子具有不同藥理活性，但是還沒有直接實驗證據(jù)揭示與本研究靶點的相關性。

表1 候選化合物及陽性對照Glide分子對接打分值與MM/GBSA結(jié)合自由能

表2 排名前4命中化合物的ADME參數(shù)預測值

圖2 虛擬篩選流程

2.2 分子動力學模擬結(jié)果分析

計算分子動力學模擬50 ns時長軌跡的均方根偏差(root mean square deviation，RMSD)，以此檢驗體系穩(wěn)定性，結(jié)果見圖3。在模擬運行35 ns時長之后，所有體系的RMSD的平均值為0.5～0.7 nm，且漲落范圍都穩(wěn)定在0.1 nm以下，均達到了平衡狀態(tài)。

圖3 靶蛋白EBOV-sGP的氨基酸骨架原子隨時間變化的RMSD值

通過計算化合物與受體蛋白復合物氨基酸骨架原子在分子動力學模擬過程中的均方根波動值(root mean square fluctuation，RMSF)，分析蛋白質(zhì)中每個氨基酸殘基的柔性，評估不同小分子抑制劑對EBOV-sGP每個氨基酸殘基波動的影響。發(fā)現(xiàn)4個抑制劑復合物體系與toremifene復合物體系都有相似的RMSF數(shù)值分布趨勢，結(jié)果見圖4。在鉸鏈區(qū)(290～310)、催化loop區(qū)(186～201、99～104)、Activationβ-sheet(514～520、62～68)這些結(jié)合位點內(nèi)的殘基RMSF都明顯較小，說明MOL006834、MOL000174、MOL002192、MOL012524一樣都可以與EBOV-sGP結(jié)合口袋形成穩(wěn)定的相互作用，穩(wěn)定蛋白結(jié)合口袋的殘基。除此之外，復合物體系(MOL0 06834、MOL000174、MOL002192、MOL012524、toremifene)與靶蛋白空載體(in an inactive or unbound state，APO form)相比較，含配體分子的蛋白二元結(jié)構(gòu)比空蛋白在活性中心關鍵氨基酸殘基520～550區(qū)間具有更小的RMSF。

圖4 配體-蛋白復合物氨基酸骨架原子的RMSF值波動圖

2.3 結(jié)合自由能計算結(jié)果分析

從35～50 ns分子動力學軌跡中提取150幀軌跡文件，計算MM/PBSA結(jié)合自由能及各種能量貢獻情況，結(jié)果見表3。范德華作用能與靜電勢能最有利于4種命中化合物與EBOV-sGP結(jié)合，表面溶劑化作用也有利于結(jié)合，但與范德華作用能和靜電勢能的數(shù)值相比作用較弱，而極性溶劑化作用不利于體系結(jié)合。表明范德華作用成分能在整個作用力占有主導地位，范德華作用能成分極大地促進了抑制劑與靶標蛋白的結(jié)合，決定了結(jié)合自由能的最終排名，這也是MOL006834、MOL000174與靶蛋白EBOV-sGP結(jié)合自由能具有優(yōu)勢效應的關鍵因素。

表3 命中化合物與靶標EBOV-sGP的MM/PBSA結(jié)合自由能

MOL006834與MOL000174來源于相同的化學結(jié)構(gòu)母核，空間結(jié)構(gòu)相似，因此范德華作用能、表面溶劑化作用、結(jié)合自由能均非常接近；但是MOL000174具有較強的極性溶劑化負效應，具體數(shù)值高達205.799 kJ/mol，非常不利于配體與受體結(jié)合，因此與蛋白質(zhì)結(jié)合過程中由于溶劑中的極性溶劑化相互作用很大程度上抵消了靜電勢能與范德華力之和，不利于導致最終的結(jié)合自由能排名次于MOL006834，位居第二。

同樣，MOL002192與MOL012524也來源于相同的化學結(jié)構(gòu)母核，空間結(jié)構(gòu)相似，因此范德華作用能、表面溶劑化作用、結(jié)合自由能均非常接近。MOL012524具有較強的極性溶劑化負效應，具體數(shù)值為156.173 kJ/mol，不利于配體與受體結(jié)合，雖然靜電勢能也是-53.094 kJ/mol，但是與蛋白質(zhì)結(jié)合過程中由于溶劑中的極性溶劑化相互作用終究抵消了靜電勢能與范德華力之和。因此MOL012524具有最弱的范德華作用正效應與較強的極性溶劑化負效應，導致最終的結(jié)合自由能不占優(yōu)勢，排名第四。

結(jié)合自由能數(shù)值越小，表示結(jié)合所需能量越低，越有利于配體與蛋白的結(jié)合。比較4種化合物與EBOV-sGP的結(jié)合自由能，可得出候選抑制劑分子結(jié)合靶標強弱關系如下：MOL006834>MOL000174>MOL002192>MOL012524，均優(yōu)于陽性對照分子toremifene。

2.4 候選化合物與靶標的相互作用位點分析

提取復合物分子動力學模擬50ns穩(wěn)定后最后一幀構(gòu)象，觀察命中化合物與EBOV-sGP的空間結(jié)合構(gòu)象。由圖5可見，靶蛋白EBOV-sGP與MOL006 834之間的相互作用不僅有范德華力、氫鍵，還有π-Anion相互作用、π-Sulfur相互作用、π-Alkyl相互作用、π-π堆積作用。MOL006834與氨基酸殘基Met548、Arg200、Pro533、Pro537、Ala539、Asn204形成常規(guī)氫鍵，并且同樣與Glu201、Asn204均形成一個碳氫鍵。與Ala101、Pro198、Pro202分別形成一個π-Alkyl相互作用，與Glu100形成π-Anion相互作用，與Met548形成π-Sulfur相互作用，與Tyr517形成π-π堆積作用。這些作用力使結(jié)合能降低，親和力增加，候選單體與EBOV-sGP蛋白的活性部位緊密結(jié)合，形成穩(wěn)定的復合物，就可以有效阻止外來病毒粒子附著并侵入宿主細胞。

圖5 靶蛋白EBOV-sGP與MOL006834相互作用

由圖6可見，靶蛋白EBOV-sGP與MOL000174之間的相互作用不僅有氫鍵、范德華力，還有π-π堆積作用、π-Alkyl相互作用。MOL000174與氨基酸殘基Arg64、Arg164、Arg200、Glu201、Thr520形成常規(guī)氫鍵，并且同樣與Gly67、Ala101、Arg200、Thr519、Glu201、Leu199形成碳氫鍵，還與Tyr517形成π-π堆積作用，與Ala101形成π-Alkyl相互作用。這些作用力使結(jié)合能降低，親和力增加，候選單體與EBOV-sGP活性空腔部位緊密結(jié)合，形成穩(wěn)定的復合物，就可以有效阻止外來病毒粒子附著并侵入宿主細胞。但是還可以觀察到氨基酸殘基Arg64與Glu100之間存在著不利相互作用，不利于構(gòu)象穩(wěn)定。

圖6 靶蛋白EBOV-sGP與MOL000174相互作用

由圖7可知，靶蛋白EBOV-sGP與MOL002192之間的主要作用力是范德華力、氫鍵、π-cation相互作用、amide-π堆積作用、π-Alkyl相互作用。MOL002192與靶標氨基酸殘基Ser195、Asp192、Arg200、Thr520、Asn61、His197均形成氫鍵，其中Asn61形成了2個常規(guī)氫鍵，Ser195分別形成了常規(guī)氫鍵和碳氫鍵各1個。與堿性氨基酸殘基Arg200、Arg64均形成π-anion相互作用，與氨基酸殘基Pro202形成π-Alkyl相互作用，并且His197形成amide- π堆積作用，與這些作用力促進復合物的形成。

圖7 靶蛋白EBOV-EBOV-sGP與MOL002192相互作用

由圖8可見，靶蛋白EBOV-sGP與MOL012524之間的相互作用不僅有范德華力、氫鍵，還有π-Alkyl相互作用。MOL012524與氨基酸殘基Glu100、Arg200、Glu201、His197形成常規(guī)氫鍵，同樣并且與Arg164、Asn98、Tyr99均形成一個碳氫鍵，與Thr519形成了一個π-Donor氫鍵，其中Glu100形成了2個常規(guī)氫鍵。與Arg64、Ala101分別形成一個π-Alkyl相互作用。這些作用力使結(jié)合能降低，親和力增加，候選單體與EBOV-sGP蛋白的活性部位緊密結(jié)合，形成穩(wěn)定的復合物，就可以有效阻止外來病毒粒子附著并侵入宿主細胞。

圖8 靶蛋白EBOV-sGP與MOL012524相互作用

氫鍵在生物大分子中起著重要作用，氫鍵是蛋白質(zhì)與配體相互作用的主要推動力之一，在穩(wěn)定蛋白質(zhì)-配體復合物方面起著不可或缺的作用[16]。Discovery Studio Visualizer生成的結(jié)果顯示，陽性對照分子toremifene與氨基酸殘基Glu100、Thr520形成氫鍵相互作用，產(chǎn)生的氫鍵數(shù)目為2個；然而4個候選單體均與EBOV-sGP 氨基酸殘基Glu201或Asn61形成氫鍵相互作用，甚至形成雙氫鍵，這些殘基均是候選單體的氫鍵受體，并且常規(guī)氫鍵總數(shù)均大于4個(其中氫鍵數(shù)MOL006834為6個，MOL0 00174為7個，MOL002192為6個，MOL012524為5個)。借助g-hbond程序統(tǒng)計氫鍵數(shù)目，計算50 ns模擬過程中EBOV-sGP與抑制劑分子間的氫鍵數(shù)目以及氫鍵距離或角度的分布。在Gromacs中判斷氫鍵要滿足兩個條件，距離和角度，g-hbond統(tǒng)計了滿足距離限制的所有原子對，并且也同時滿足角度限制，即可判定為氫鍵，統(tǒng)計結(jié)果如圖9所示。總體而言，與陽性對照分子toremifene相比，命中的候選分子與EBOV-sGP結(jié)合氫鍵數(shù)目更加具有優(yōu)勢，均具有更強的相互作用，推測可能具有更高的抑制作用。

圖10展示了陽性對照分子toremifene與靶標的結(jié)合構(gòu)象，可以觀察到toremifene與EBOV-sGP 氨基酸殘基Val66、Glu515均形成π-sigma相互作用，氨基酸殘基Ala101同時與配體分子的兩個苯環(huán)形成π-alkyl相互作用，與堿性氨基酸Arg64形成π-Cation堆積作用，3個Leu氨基酸殘基184、186、558同時與配體分子末端氧原子形成alkyl相互作用，這些作用力促進復合物的形成。同樣在4個候選單體與靶蛋白的相互作用分析中，MOL006834、MOL002192均與氨基酸殘基Pro202形成π-Alkyl相互作用，MOL000174、MOL012524與堿性氨基酸Ala101也同樣均形成π-Alkyl相互作用，這些相互作用在候選配體分子與EBOV-sGP單體空腔的殘基結(jié)合中也起了輔助有利作用。

圖10 靶蛋白EBOV-sGP與toremifene相互作用

3 討論與結(jié)論

Ahmad課題組于2017年首次基于EBOV-GP跨膜糖蛋白三聚體空腔的三維結(jié)構(gòu)為靶標，從Mcule化合物庫(https://mcule.com/database/)進行大規(guī)模虛擬篩選抑制劑，采用Auto Dock Vina與Flex-X兩種分子對接軟件，并結(jié)合網(wǎng)上在線分析工具admetSAR(http://www.admetexp.org)進行ADMET過濾，找到了三個命中分子[17]。從2000年以來，不斷有埃博拉病毒疫情暴發(fā)，抗病毒藥效分子的研發(fā)仍是全球醫(yī)藥研究熱點。因此本研究為了從中藥天然產(chǎn)物庫中挖掘抗病毒活性分子，借鑒了此文章的研究思路，靶標從EBOV-GP三聚體改為單體，并重新設計了多輪虛擬篩選的實驗方法，篩選策略使用成本更高的算法，盡可能提高預測的精度，降低虛擬篩選假陽性率，以及更好地預測結(jié)合構(gòu)象。除此之外，分子動力學模擬和結(jié)合自由能計算也是研究蛋白質(zhì)與抑制劑相互作用強度的重要工具，采用GROMACS軟件進行分子動力學模擬，并使用 MM/PBSA算法預測體系穩(wěn)定后的結(jié)合自由能，這類方法不僅能在原子和分子的層次上給出蛋白質(zhì)和抑制劑的相互作用細節(jié)，而且還能便于闡明蛋白質(zhì)和抑制劑復合體的結(jié)構(gòu)親和能關系。

為了有效干擾EBOV-sGP構(gòu)象轉(zhuǎn)變，本研究基于EBOV-GP單體活性空腔為靶標，從中藥天然產(chǎn)物庫中挖掘抑制劑分子，有助于阻斷宿主-病毒膜融合。首先遵循Lipinski五原則、Verber三原則，過濾掉違反類藥性原則的化合物；然后采用Glide軟件進行三種不同精度的分子對接計算，即高通量HTVS、標準精度SP、超精密XP等三種模式，逐步縮小候選化合物庫的規(guī)模；隨后采用Schrodinger的QikProp程序預測命中分子ADMET性質(zhì)，較admetSAR進行過濾策略，新增了肝毒性HERG_IC50、口服生物利用度等兩項藥代動力學與藥物安全性分析的關鍵指標，并結(jié)合Lipinski Ro5、Verber Ro3預測值，做出相對全面的類藥性評估；接著在分析靶標蛋白與命中分子結(jié)合方式過程中，使用高性能GPU加速技術模擬在真實生理鹽水溶劑條件下生物大分子的運動狀態(tài)，進行了50 ns時長的分子動力學模擬，并使用了MM/GBSA與MM/PBSA兩種算法預測結(jié)合自由能，提升富集精準度。最后根據(jù)XP-Glide 評分、RMSD、RMSF、MM/GBSA 結(jié)合能、MM/PBSA結(jié)合能、ADME預測，命中分子MOL006834、MOL000174、MOL002192、MOL012524表現(xiàn)出對 EBOV-sGP活性位點顯示出極好的結(jié)合親和力和良好的類藥性；從原子與分子層次上預測命中分子MOL006834、MOL000174、MOL002192、MOL012524與跨膜糖蛋白的真實結(jié)合構(gòu)象，包括氫鍵、疏水相互作用、離子鍵、π-π堆積力等，與關鍵氨基酸殘基顯示出極好的相互作用。這四個命中分子均可在商業(yè)上購買獲得，也可作為EBOV-sGP抑制劑先導化合物進行開發(fā)利用，以及藥物活性實驗驗證。需要指出的是，埃博拉病毒具有高致死率和高傳播性，相關生物安全要求別極其嚴格，抗病毒活性測試實驗需要在生物安全4級實驗室內(nèi)進行,國內(nèi)具備此級別的實驗室并不多，制約了有關埃博拉病毒的基礎研究。因此本研究成果僅局限于藥物開發(fā)的早期階段，缺乏抗病毒活性驗證、ADMET評估等實驗數(shù)據(jù)支撐，后期需要聯(lián)合相關的病毒學研究團隊，借助專業(yè)的病毒學研究平臺開展安全有效的生物學驗證。