史海龍,程 怡,黃 月,馮雪松,王月雯,黃 峰,晁 旭
基于計算機模擬技術從中藥天然產(chǎn)物庫中挖掘色氨酸羥化酶1抑制劑
史海龍,程 怡,黃 月,馮雪松,王月雯,黃 峰,晁 旭*
陜西中醫(yī)藥大學,陜西 西安 712046
基于計算機模擬技術從中藥天然產(chǎn)物庫中挖掘色氨酸羥化酶1(tryptophan hydroxylase-1,TPH1)抑制劑。搭建分子對接、類藥性篩選、藥動學預測、分子動力學模擬于一體的藥物篩選平臺,從中藥系統(tǒng)藥理學數(shù)據(jù)庫與分析平臺(traditional chinese medicine systems pharmacology database and analysis platform,TCMSP)數(shù)據(jù)庫中挖掘TPH1抑制劑。川貝酮堿(TCMSP_ID:MOL009572)實驗數(shù)據(jù)均表現(xiàn)出良好的TPH1抑制活性及類藥性,預測命中分子可有效抑制胃腸道相關TPH1活性,但對中樞神經(jīng)系統(tǒng)相關TPH2活性抑制作用較小;通過酶活抑制實驗驗證抑制效應,并借助分子動力學模擬動態(tài)解析結合能變化及各項能量分布。整合多種虛擬篩選技術挖掘到緩解腸易激綜合征相關痛癥的TPH1抑制劑-川貝酮堿MOL009572。
分子對接;分子動力學模擬;色氨酸羥化酶;抑制劑;川貝酮堿
5-羥色氨酸(5-hydroxytryptophan,5-HT)信號通路在內(nèi)臟痛覺高敏感、腸動力異常等效應產(chǎn)生中扮演著重要角色[1],科學家們期待以此為突破口,尋求針對腸易激綜合征、潰瘍性結腸炎等胃腸道疾病的治療策略。色氨酸羥化酶(tryptophan hydroxylase,TPH)可以使色氨酸的羥基化,使色氨酸轉化為5-HT。這是5-HT生物合成過程中的起始步及限速步,因此TPH是治療或緩解5-HT相關腸易激綜合征(irritable bowel syndrome,IBS)痛癥藥物的關鍵靶標[2]。哺乳動物TPH普遍存在2種亞型,其中人類TPH1與TPH2基因分別位于11、12號染色體上,序列相似性高達71%[3]。TPH1在胃腸道等消化道上皮組織中表達,TPH2主要在中樞神經(jīng)系統(tǒng)的神經(jīng)細胞中表達[4],因此利用TPH1表達的組織特異性研發(fā)抑制劑,阻斷消化系統(tǒng)內(nèi)5-HT合成,可發(fā)揮IBS鎮(zhèn)痛止痛的療效。目前已有TPH1抑制劑類藥物上市,但臨床用藥中仍然存在諸多問題,如長期服用可能導致毒副作用甚至耐藥性,或引發(fā)腦部中樞5-HT合成阻斷[5],因此需要研發(fā)新型TPH1抑制劑。
近幾年研究相繼解析出高分辨率TPH蛋白晶體結構,推動了基于靶標活性口袋結構的計算機輔助藥物分子設計,然而國內(nèi)外IBS適應癥藥物研發(fā)主要圍繞著TPH1為靶標開發(fā)新型單一抑制劑/拮抗劑,未涉及到TPH2酶活抑制的負效應。因此本研究搭建類藥性篩選、藥動學參數(shù)預測、雙靶標分子對接、分子動力學模擬于一體的藥物篩選平臺,從中藥天然產(chǎn)物庫中挖掘TPH1抑制劑。
分子對接均采用Schrodinger 2018的Glide模塊、ADME參數(shù)計算采用Schrodinger 2018的QikProp模塊,隱性溶劑環(huán)境的結合自由能計算采用分子力學/廣義波恩表面積模型MM/GBSA[6],顯性溶劑環(huán)境的結合自由能計算采用分子力場/泊松玻爾茲曼表面積模型MM/PBSA[7],分子動力學模擬采用Gromacs 2019與Gaussian 09程序包。
人源性的重組表達TPH1、TPH2由陜西中醫(yī)藥大學針藥結合省級重點實驗室朱先偉博士提供;LX-1031(CAS 945976-76-1,批號637958)購于abcam中國上海分公司;川貝酮堿(CAS 103530-47-8,批號20190718)購于上海同田生物;ELX808酶標儀(美國伯騰公司)。
首先從PDB數(shù)據(jù)庫下載的靶蛋白晶體結構,TPH1選擇3HF8(PDB_ID),TPH2選擇4V06(PDB_ID),見圖1;然后采用Schrodinger的Protein Preparation Wizard模塊預處理受體蛋白結構;隨后采用Protein Grid Generation模塊準備靶標TPH1/TPH2蛋白活性口袋并產(chǎn)生格點文件,原結晶配體坐標作為格點盒子中心。由于篇幅限制,僅展示靶標TPH1與原共結晶配體LP533401結合模式及關鍵殘基,見圖2。

圖1 靶標受體蛋白TPH1和TPH2活性口袋

圖2 晶體復合物3HF8中TPH1與原配體LP-533401結合模式
中藥天然產(chǎn)物數(shù)據(jù)庫TCMSP[8](https:// old.tcmsp-e.com/index.php)13 144個化合物作為配體庫。采用Schrodinger的LigPrep模塊對各分子進行加氫、能量最小化、幾何優(yōu)化等處理。在力場OPLS3下,計算pH為7.0±2.0的化合物離子狀態(tài),搜索合理構象作為對接起始構象,形成可能的互變異構體。由于LX-1031有效抑制TPH1酶活性,因此可作為陽性對照,用于5-HT過表達相關疾病的藥物研發(fā),如腹瀉型IBS[9]。以同樣方式預處理LP533401、LX-1031,作為本實驗對照組。
基于類藥性規(guī)則“Lipinski Ro5[10]”和“Verber Ro3[11]”對中藥分子配體庫進行首輪篩選,挑選出8987個符合類藥性規(guī)則的命中分子。
基于受體結構的雙靶標TPH1/TPH2分子對接進行第2輪篩選。首先針對靶標TPH1對上一步篩選的化合物庫進行高通量虛擬篩選(high throughput virtual screening docking,HTVS- docking),并根據(jù)Glide Score對HTVS-docking的輸出化合物進行排名。然后提取最高得分前35%化合物,進行靶標TPH1標準精度(standard precision docking,SP-docking)模式分子對接,根據(jù)Glide Score對SP-docking的輸出化合物再次進行排名。隨后提取最高得分前40%化合物,進行TPH1/TPH2雙靶標高精度(extra precision docking,XP-docking)模式分子對接,以TPH1的Docking score≤?33.49 kJ/mol為篩選閾值,選取與TPH1有較低結合自由能的分子,再次選取與TPH2有較高Docking score的化合物,二者交集最終作為候選分子。
基于ADME藥動學參數(shù),采用Schrodinger的QikProp模塊進行第3輪篩選,從相對分子質量(molecular weight,W)、氫鍵供體(hydrogen bond donor,HB donor)、氫鍵配體(hydrogen bond acceptor,HB acceptor)、溶劑可及表面積(solvent accessible surface area,SASA)、K+通道阻斷相關log IC50系數(shù)(log IC50for blockage of HERG K+channels,logHERG)、水溶性系數(shù)(aqueous solubility,logS)、脂水分配系數(shù)[octanol/water partition coefficient,log(O/W)]、口服生物利用度(oral absorption,OB)、化合物的極性表面積(PSA,Van der Waals surface area of polar nitrogen and oxygen atoms and carbonyl carbon atoms)等評估,篩選具有良好類藥性質的化合物。具體篩選規(guī)則為W可接受范圍小于650、HB donor數(shù)值可接受范圍不大于5、HB acceptor數(shù)值可接受范圍不大于10、SASA數(shù)值可接受范圍為300~1000、K+通道阻斷相關log IC50(數(shù)值推薦區(qū)間低于-5)、logS數(shù)值可接受范圍為?7.0~0.5)、log(O/W)數(shù)值可接受范圍為?2~6.5、OB低于50代表較差,高于80代表較好,PSA數(shù)值可接受范圍為7.0~200.0。
基于Prime-MM/GBSA結合自由能預測進行第4輪篩選。由于高精度分子對接函數(shù)計算仍然無法估算配體-受體復合物結合的溶劑化效應因素,因此使用隱性溶劑模型下MM/GBSA算法優(yōu)化排序,進一步提高結合自由能的計算精度。
為了從動力學和熱力學角度動態(tài)分析靶標TPH1與命中分子互作模式,使用GROMACS進行分子動力學模擬。蛋白拓撲參數(shù)取自Amber的FF03力場,由PDB2GMX工具獲得;抑制劑的各原子電荷先用Gaussian 09在B3LYP/6-311G**基組水平計算靜電勢,再用Amber中的RESP電荷擬合程序算出。將復合物放入1 nm×1 nm×1 nm的立方體TIP3P水盒中,并加入Na+和Cl?,中和該系統(tǒng)到0.15 mol/L NaCl溶液濃度,隨后進行能量最小化和預平衡,最后在26.85 ℃(300 K)進行50 ns的分子動力學模擬,時間間隔為2 fs,每10 ps保存一次能量和坐標文件。從25~50 ns穩(wěn)定模擬軌跡中提取250幀軌跡文件,計算MM/PBSA結合自由能,并且提取TPH1- MOL009572、TPH1-LX-1031兩個體系分子動力學模擬50 ns最后一幀構象,解析受體-配體結合模式。
命中分子MOL009572進行TPH酶活抑制的驗證性實驗。酶活測定緩沖液包含:60 μmol/L色氨酸溶液、300 mol/L 6-甲基-四氫蝶呤、200 mmol/L硫酸銨溶液、7 mmol/L DTT、25 μg/mL過氧化氫酶、25 μmol/L硫酸亞鐵銨、50 mmol/L MES緩沖液(pH 7.0)。25 nmol/L人源性TPH1/TPH2被添加到不同濃度梯度的MOL009572溶液中開始反應。使用4 mm的石英比色皿熒光分光光度法定量測定酶活性,激發(fā)光波長為300 nm,發(fā)射光波長為330 nm,分別于狹縫寬度在3 nm處激發(fā)和5 nm處發(fā)射。
TCMSP數(shù)據(jù)庫中13 144個中藥單體化合物進行多輪虛擬篩選,篩選流程見圖3,挑選出排名前10的命中分子,它們均具備有較強的TPH1結合自由能、較弱的TPH2結合自由能、良好的藥動學參數(shù)、較高的OB等,如表1、2所示。進一步參照LX-1031的Prime-MM/GBSA結合自由能,以?57.73 kJ/mol作為閾值,篩選得到2個抑制劑分子、MOL009572(川貝酮堿,chuanbeinone)、MOL009586(異浙貝母堿,isoverticine),二者均是來源于川貝母的等甾體生物堿,結構式見圖4。
計算分子動力學模擬50 ns時長軌跡的RMSD值,以此檢驗體系穩(wěn)定性,結果見圖5。在模擬運行15 ns之后,所有體系的RMSD值均無大幅波動,均達到了平衡狀態(tài)。從25~50 ns分子動力學軌跡中提取250 幀軌跡文件,計算MM/PBSA結合自由能及各種能量貢獻情況,結果見表3。結合自由能強弱關系為MOL009572>LX-1031>MOL009586>LP533401,其中范德華力最有利于命中分子與靶標結合,并起主導作用;靜電勢能、表面溶劑化作用也有利于結合,但作用力較弱,而溶劑中極性溶劑化作用不利于體系結合。在LP-533401與靶標結合中,靜電勢能、范德華力均有利于配體和受體的結合,靜電勢能處于主導地位;但是極性溶劑化作用力較強,不利于體系結合;在候選分子MOL009586中,與靶標的范德華力較弱,而溶劑中的極性相互作用較高,2個因素均不利于體系結合,造成與TPH1結合自由能數(shù)值偏高。基于以上分析,選擇MOL009572作為TPH1的抑制分子。

圖3 虛擬篩選流程及結果

表1 命中分子(排名前10) 雙靶標TPH1/TPH2對接結合能及MM/GBSA結合自由能
計算15~50 ns TPH1蛋白骨架原子的均方根波動值(root mean square fluctuation,RMSF),結果見圖6。體系TPH1-MOL009572、TPH1-LX-1031中蛋白骨架原子波動變化趨勢相似,氨基酸殘基222~225區(qū)段和361~371區(qū)段出現(xiàn)較大波動,無法形成穩(wěn)定的成鍵相互作用,不利于受體-配體復合物的結合,但是2個體系的整體性波動幅度均小于空載體蛋白靶標,說明MOL009572、LX-1031與靶標的結合,有利于靶標蛋白空間結構穩(wěn)定。除此之外,還觀察到2個體系在氨基酸殘基270~300區(qū)段、310~320區(qū)段出現(xiàn)了較小波動(RMSF<0.1 nm),TPH1活性口袋周圍氨基酸殘基(Leu236、Pro238、Arg257、Pro262、Phe263、Tyr264、Thr265、Pro266,Glu267、Pro268、Asp269、Cys271、Thr310、Ile366、Thr367、Thr368,見表4)RMSD波動值均小于空受體蛋白APO、LX-1031相應的殘基位點,推測MOL009572與TPH1之間形成了氫鍵、疏水相互作用、π-π堆積作用等,將MOL009572包裹在蛋白活性口袋內(nèi),形成更穩(wěn)定的受體-配體二元復合物。

表2 命中分子(排名前10)的ADME參數(shù)預測值

圖4 命中分子和陽性對照分子的化學結構
如圖7所示,分別為MOL009572對酶TPH1、TPH2活性抑制實驗和LX-1031對酶TPH1、TPH2活性抑制實驗,這4組實驗數(shù)據(jù)均由酶動力學實驗方法測量得到。圖7-a的實驗數(shù)據(jù)計算得出,MOL009572與TPH1、LX-1031與TPH1二者的體外IC50相同,均為210 μmol/L。圖7-b的實驗數(shù)據(jù)計算得出,LX-1031與TPH2的體外IC50160 μmol/L,MOL009572與TPH2未觀察到抑制效果。需要提出說明的是,以上每個濃度均做3次平行重復,分別取酶抑制活性平均值。此外從圖7-c實驗數(shù)據(jù)作Line-weaver-Burk曲線,結果表明MOL009572顯示出競爭性抑制的特點。

圖5 分子動力學模擬過程中復合物的TPH1氨基酸骨架原子隨時間變化的RMSD值

表3 MM/PBSA結合自由能

圖6 分子動力學模擬過程中復合物的TPH1氨基酸骨架原子隨時間變化的RMSF值
Discovery Studio Visualizer分析MOL009572與TPH1互作的主要作用力是氫鍵、π-Alkyl堆積作用,見圖8。MOL009572與氨基酸殘基Tyr264形成常規(guī)氫鍵,與Thr265形成碳氫鍵,并且與在Val232、Leu236、Pro268形成Alkyl相互作用,與Phe241、Tyr235、Phe318、His272形成π-Alkyl堆積作用,這些作用力促進受體-配體相互結合。如圖9所示,TPH1與LX-1031的相互作用不僅有氫鍵,還有π-π堆積作用。LX-1031與殘基Glu317、Ser336、Arg257、Thr265形成常規(guī)氫鍵,與殘基Tyr235、Arg121形成碳氫鍵,同時還與Tyr235、Phe313、His272形成π-π堆積作用,與Ile366、Pro268形成π-Alkyl或Alkyl相互作用。這些作用力使結合能降低,LX-1031親和力增加,與TPH1活性部位緊密結合,形成穩(wěn)定的受體-抑制劑復合物,降低TPH1酶的催化活性,從而有效抑制5-HT合成。

表4 分子動力學模擬過程中TPH1活性口袋部位關鍵殘基RMSF值

a-對人源性TPH1的抑制效果 b-對人源性TPH2的抑制效果 c-不同濃度抑制劑MOL009572的初速度的Lineweaver-Burk曲線
借助g-hbond程序計算TPH1與MOL009572、LX-1031之間的氫鍵數(shù)量隨時間的變化,統(tǒng)計結果如圖10所示。總體而言,與LX-1031相比,MOL009572與TPH1結合氫鍵數(shù)量近似;其中MOL009572氫鍵數(shù)量峰值為3,均值為3;LX-1031氫鍵數(shù)量峰值為4,均值為2,說明氫鍵促使復合物形成穩(wěn)定構象。在分子動力學模擬45 ns之后,LX-1031與TPH1至少存在3次氫鍵相互作用的接觸,而LX-1031與TPH1僅存在至少2次氫鍵相互作用的接觸,提示MOL009572比LX-1031與靶標TPH1結合更具優(yōu)勢。
四川大學華西醫(yī)學院何菱研究團隊等使用分子動力學模擬、結合自由能預測、能量項拆分及丙氨酸掃描技術,研究TPH1抑制劑苯丙氨酸衍生物的抑制作用機制[12]。結果表明,萘取代基的同分異構體可以與靶標形成不同的結合模式,但是導致相同酶活抑制效果。隨后該團隊使用基于配體結構的闡明一系列已知TPH1抑制劑的定量構效關系,并提出多元藥效團模型的研究策略[13],即整合多種TPH1抑制劑復合物晶體結構構建藥效團模型。該模型成功預測了32種苯丙氨酸取代基衍生物的酶活抑制常數(shù)。目前國內(nèi)外內(nèi)臟痛覺高敏感的藥物研發(fā)策略大體類似于上述案例,主要考慮TPH1單一靶標的抑制劑,未涉及到TPH2酶活抑制引起的負效應,篩選策略均較為單一,預測精度有待提升。

圖8 TPH1與MOL009572相互作用

圖9 TPH1與LX-1031相互作用

圖10 分子動力學模擬過程中配體小分子與受體蛋白活性口袋形成的氫鍵數(shù)目
基于上述考慮,本研究提高從挖掘傳統(tǒng)中藥庫藥效成分為突破口,制定了一套整合類藥性篩選、藥動學參數(shù)預測、雙靶標分子對接、結合自由計算的研究策略;對TCMSP數(shù)據(jù)庫進行系統(tǒng)篩選,找到一種來源于中藥川貝母的川貝酮堿(chuanbeinone,TCMSP-ID:MOL009586),隨后進行酶活抑制實驗驗證,初步證明其可抑制胃腸道相關TPH1酶活性,并降低了中樞神經(jīng)系統(tǒng)相關TPH2酶活抑制的負效應。最后采用分子動力學模擬技術從原子水平上深度剖析命中分子與靶標TPH1相互作用方式,并且分別在隱性溶劑環(huán)境下和顯性溶劑環(huán)境下計算結合自由能。本研究挖掘IBS相關痛癥的中藥藥效分子,可為研發(fā)新型TPH1抑制劑提供參考。然而,本研究還處于藥物開發(fā)的前期階段,尚缺乏藥動學評估等實驗數(shù)據(jù)支撐,后期將聯(lián)合實驗藥理學團隊,開展相關生物學活性驗證及安全性評價。
利益沖突 所有作者均聲明不存在利益沖突
[1] Okaty B W, Commons K G, Dymecki S M.Embracing diversity in the 5-HT neuronal system [J]., 2019, 20(7): 397-424.
[2] Coates M D, Tekin I, Vrana K E,.Review article: The many potential roles of intestinal serotonin (5-hydroxytryptamine, 5-HT) signalling in inflammatory bowel disease [J]., 2017, 46(6): 569-580.
[3] Walther D J, Bader M.A unique central tryptophan hydroxylase isoform [J]., 2003, 66(9): 1673-1680.
[4] del Colle A, Israelyan N, Gross Margolis K.Novel aspects of enteric serotonergic signaling in health and brain-gut disease [J]., 2020, 318(1): G130-G143.
[5] Scotton W J, Hill L J, Williams A C,.Serotonin syndrome: Pathophysiology, clinical features, management, and potential future directions [J]., 2019, 12: 1178646919873925.
[6] Zhang X H, Perez-Sanchez H, Lightstone F C.A comprehensive docking and MM/GBSA rescoring study of ligand recognition upon binding antithrombin [J]., 2017, 17(14): 1631-1639.
[7] Kongsted J, Ryde U.An improved method to predict the entropy term with the MM/PBSA approach [J]., 2009, 23(2): 63-71.
[8] Ru J L, Li P, Wang J N,.TCMSP: A database of systems pharmacology for drug discovery from herbal medicines [J]., 2014, 6: 13.
[9] Camilleri M.LX-1031, a tryptophan 5-hydroxylase inhibitor, and its potential in chronic diarrhea associated with increased serotonin [J]., 2011, 23(3): 193-200.
[10] Lipinski C A, Lombardo F, Dominy B W,.Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings [J]., 2001, 46(1/2/3): 3-26.
[11] Veber D F, Johnson S R, Cheng H Y,.Molecular properties that influence the oral bioavailability of drug candidates [J]., 2002, 45(12): 2615-2623.
[12] Ouyang L, He G, Huang W,.Combined structure-based pharmacophore and 3D-QSAR studies on phenylalanine series compounds as TPH1inhibitors [J]., 2012, 13(5): 5348-5363.
[13] Zhong H, Huang W, He G,.Molecular dynamics simulation of tryptophan hydroxylase-1: Binding modes and free energy analysis to phenylalanine derivative inhibitors [J]., 2013, 14(5): 9947-9962.
Exploring potential TPH1 inhibitors from traditional Chinese medicine natural product database based on computer simulation technology
SHI Hai-long, CHENG Yi, HUANG Yue, FENG Xue-song, WANG Yue-wen, HUANG Feng, CHAO Xu
Shaanxi University of Chinese Medicine, Xi’an 712046, China
To discover the active molecules with inhibitory activity to tryptophan hydroxylase-1 (TPH1) from traditional Chinese medicine natural product database based on computer simulation technology.A drug screening platform, that integrates molecular docking, drug-like screening, pharmacokinetic prediction, and molecular dynamics simulation, was built to discovery TPH1 inhibitors from the TCMSP database.Chuanbeinone (TCMSP_ID: MOL009572) showed good inhibitory activity to TPH1 and drug-likeness, which can effectively inhibit gastrointestinal related TPH1 enzymes, but had weak inhibitory effect on central nervous system-related TPH2 enzyme activity.Through preliminary verification of enzyme activity inhibition, molecular dynamics simulation is used to accurately predict the free binding energy of THP1 and MOL009572 and various energy contributions.Chuanbeidone, MOL009572, an active molecule of traditional Chinese medicine that relieves irritable bowel syndrome-related pain, has been screened out through integrated multiple virtual screening technologies.
molecular docking; molecular dynamics simulation; tryptophan hydroxylase; inhibitor; chuanbeinone
R284.1
A
0253 - 2670(2022)10 - 2968 - 09
10.7501/j.issn.0253-2670.2022.10.005
2021-11-21
陜西省中醫(yī)管理局中醫(yī)藥科研課題(JCMS003);陜西中醫(yī)藥大學創(chuàng)新團隊項目(2019-Y505)
史海龍,男,博士,美國UNMC博士后,副教授,研究方向為計算機輔助藥物設計及生物大分子模擬。E-mail: shl112@sntcm.edu.cn
通信作者:晁 旭,博士,教授,研究方向為中醫(yī)藥防治消化系統(tǒng)腫瘤。
[責任編輯 王文倩]