林 峰 付新梅 王 超 蔣思宇 王景輝張述偉 楊 凌 李 燕,*
(1大連理工大學(xué)化工學(xué)院,遼寧大連116024;2大連理工大學(xué)精細(xì)化工國家重點(diǎn)實(shí)驗(yàn)室,遼寧大連116024;3中國科學(xué)院大連化學(xué)物理研究所,藥物資源開發(fā)研究組,遼寧大連116023)
3C-like蛋白酶抑制劑的構(gòu)效關(guān)系、分子對(duì)接和分子動(dòng)力學(xué)
林 峰1付新梅2,*王 超1蔣思宇1王景輝1張述偉1楊 凌3李 燕1,*
(1大連理工大學(xué)化工學(xué)院,遼寧大連116024;2大連理工大學(xué)精細(xì)化工國家重點(diǎn)實(shí)驗(yàn)室,遼寧大連116024;3中國科學(xué)院大連化學(xué)物理研究所,藥物資源開發(fā)研究組,遼寧大連116023)
3C-like蛋白酶是中東呼吸綜合征冠狀病毒(MERS-CoV)等其它冠狀病毒的繁殖過程中極為重要的蛋白酶。它已成為人類在抗冠狀病毒領(lǐng)域中的研究熱點(diǎn)。本文基于計(jì)算生物學(xué)方法對(duì)與MERS-CoV同屬的蝙蝠冠狀病毒HKU4(HKU4-CoV)的43個(gè)肽類3C-like蛋白酶抑制劑分子,建立三維定量構(gòu)效關(guān)系(3D-QSAR)模型。在基于配體疊合的基礎(chǔ)上,發(fā)現(xiàn)比較分子相似性指數(shù)分析法(CoMSIA)中的四個(gè)場(chǎng)組合(位阻場(chǎng)、靜電場(chǎng)、氫鍵供體場(chǎng)與氫鍵受體場(chǎng))為最優(yōu)的模型(Q2=0.522,R2ncv=0.996,R2pre=0.904;Q2:交叉驗(yàn)證相關(guān)系數(shù),R2ncv:非交叉驗(yàn)證相關(guān)系數(shù),R2pre:驗(yàn)證集分子的預(yù)測(cè)值相關(guān)系數(shù)),并借助該模型通過分子對(duì)接(docking)與分子動(dòng)力學(xué)(MD)方法闡明了配受體結(jié)合作用。實(shí)驗(yàn)結(jié)果表明:(1)基于最優(yōu)的CoMSIA模型基礎(chǔ)上的三維等勢(shì)圖形象地說明了分子基團(tuán)的位阻作用、靜電作用、氫鍵供體與氫鍵受體作用對(duì)分子生物活性的影響;(2)分子對(duì)接研究結(jié)果顯示了疏水性以及結(jié)晶水、氨基酸His166和Glu169在配體和受體結(jié)合過程中產(chǎn)生重要作用;(3)分子動(dòng)力學(xué)模擬進(jìn)一步驗(yàn)證了分子對(duì)接模型的可靠性,并發(fā)現(xiàn)了兩個(gè)新的關(guān)鍵氨基酸Ser24與Gln192,它們與配體產(chǎn)生了兩個(gè)較強(qiáng)的氫鍵。此外,根據(jù)這些結(jié)果,一些新的具有潛在抑制活性的肽類化合物作為3C-like蛋白酶抑制劑被獲得。以上結(jié)果能夠幫助深入了解3C-like蛋白酶與肽類抑制劑的作用機(jī)理,并且能夠?yàn)榻窈蟮目筂ERS-CoV藥物設(shè)計(jì)提供有價(jià)值的參考。
中東呼吸綜合征冠狀病毒;3C-like蛋白酶;肽類抑制劑;三維定量構(gòu)效關(guān)系;分子對(duì)接;分子動(dòng)力學(xué)
冠狀病毒(CoV)是一類分布廣泛的,具有重大威脅的病原體1,2。該種病毒能引起人及哺乳類動(dòng)物疾病,且被感染的人或動(dòng)物可能會(huì)成為呼吸道、腸道、肝臟、神經(jīng)系統(tǒng)疾病的攜帶者。最近出現(xiàn)的中東呼吸綜合征冠狀病毒(MERS-CoV)就是其中能感染人類并傳播的冠狀病毒之一3。該病毒于2012年9月在中東地區(qū)首次被鑒定出來,隨后逐漸擴(kuò)散,在很多國家都報(bào)道了感染病例。MERS-CoV能導(dǎo)致非常嚴(yán)重的呼吸道病癥,致死率高。截至2015年1月30日,共有956例確診病例,死亡351例,死亡率高達(dá)36.7%,遠(yuǎn)高于嚴(yán)重急性呼吸綜合征冠狀病毒(SARS-CoV),引發(fā)了全球關(guān)于MERS大流行潛在可能性的恐慌4。因此,MERS-CoV已經(jīng)成為一種全球范圍內(nèi)嚴(yán)重威脅人類生命安全的新型冠狀病毒。
冠狀病毒被國際病毒分類委員會(huì)(ICTV)分為四個(gè)屬,即α、β、γ和δ冠狀病毒屬。其中,β冠狀病毒屬有四個(gè)亞群(A-D),通過將MERS-CoV和其他β冠狀病毒基因序列對(duì)比發(fā)現(xiàn),MERS-CoV與SARS-CoV中較保守的核酸復(fù)制酶序列只有小于50%的氨基酸相同,而與從蝙蝠體內(nèi)分離到的蝙蝠冠狀病毒HKU4(HKU4-CoV)和HKU5 (HKU5-CoV)親緣關(guān)系很近,相似度分別約為75%和76.7%5-7。因此,ICTV在β冠狀病毒屬C亞群中增加MERS-CoV為一個(gè)新種,它也是該亞群中第一種能夠感染人類的β冠狀病毒。研究表明,MERS-CoV不僅可存在于蝙蝠體內(nèi),其也被發(fā)現(xiàn)能生存在卡塔爾單峰駱駝體內(nèi),而且它們的輔助蛋白在體外能夠抑制人類抗病毒的信號(hào)路徑3,8-10。除此之外,Wang等11通過研究認(rèn)為MERS-CoV的起源是HKU4-CoV或HKU5-CoV的變異,其中更多研究表明,其更傾向是來自HKU4-CoV的變異,表明從蝙蝠到人類的一種人畜共患的轉(zhuǎn)變。
3C-like蛋白酶(3CLpro)是病毒繁殖的蛋白水解過程中,兩個(gè)非常重要的病毒半胱氨酸蛋白酶之一12-14。3C-like蛋白酶在擁有典型的Leu-Gln (Ser、Ala、Gly)氨基酸序列的11個(gè)保留位點(diǎn)上將復(fù)制酶多聚蛋白斷裂,其中它們的P1位置與P2位置分別擁有完整的Gln氨基酸殘基和脂肪族氨基酸殘基15,16。迄今為止,沒有疫苗和抗病毒藥物能夠真正有效預(yù)防冠狀病毒感染人類或治療已患病的人類。因此,3C-like蛋白酶對(duì)于冠狀病毒生命周期起著至關(guān)重要的作用,使其成為一個(gè)抗MERSCoV藥物的發(fā)展很有吸引力的目標(biāo),基于3C-like蛋白酶抑制劑的研究具有廣闊的前景17-19。
隨著計(jì)算機(jī)科技的高速發(fā)展,三維定量構(gòu)效關(guān)系(3D-QSAR)作為一個(gè)十分有效、經(jīng)濟(jì)的方法,已被廣泛應(yīng)用于各個(gè)領(lǐng)域的結(jié)構(gòu)與活性、結(jié)構(gòu)與性質(zhì)關(guān)系的研究,來幫助設(shè)計(jì)新的藥物20,21。3DQSAR包括比較分子力場(chǎng)分析(CoMFA)和在其基礎(chǔ)之上發(fā)展起來的較新穎的比較分子相似性指數(shù)分析(CoMSIA)等眾多方法。這些方法認(rèn)為,配體與受體之間的相互作用取決于化合物周圍分子力場(chǎng)的差異,以定量化的分子力場(chǎng)參數(shù)作為變量,對(duì)藥物活性進(jìn)行回歸分析便可反映出藥物小分子與生物大分子之間的作用模式,進(jìn)而有選擇地設(shè)計(jì)新藥。
目前MERS-CoV的3C-like蛋白酶抑制劑還處于研究階段,但已有文獻(xiàn)證明與MERS-CoV同屬的SARS-CoV的3C-like蛋白酶抑制劑也同樣能良好地抑制MERS-CoV的繁殖22-25。本文通過上述計(jì)算方法,對(duì)新合成的與MERS-CoV更接近的HKU4-CoV的43個(gè)肽類3C-like蛋白酶抑制劑分子進(jìn)行了研究,建立了3D-QSAR模型,并探討了此類抑制劑周圍的位阻場(chǎng)、靜電場(chǎng)、氫鍵供體場(chǎng)和氫鍵受體場(chǎng)對(duì)生物活性的影響。此外,我們還利用分子對(duì)接方法來驗(yàn)證3D-QSAR模型的準(zhǔn)確性和穩(wěn)定性,之后采用分子動(dòng)力學(xué)方法對(duì)分子對(duì)接的結(jié)果進(jìn)行補(bǔ)充。它們能幫助我們更好地預(yù)測(cè)抑制劑小分子與靶點(diǎn)蛋白大分子之間可能的結(jié)合構(gòu)象,進(jìn)而預(yù)測(cè)小分子與蛋白的結(jié)合力和其生物活性。我們希望獲得的分子構(gòu)效關(guān)系以及配受體結(jié)合作用的結(jié)論能夠延伸到為今后抗MERS-CoV藥物(3C-like蛋白酶抑制劑)的設(shè)計(jì)提供理論指導(dǎo)。
2.1數(shù)據(jù)庫及生理活性
本文模型的建立及分析采用的是Tripos公司的SYBYL 6.9設(shè)計(jì)軟件包。文中所研究具有抑制3C-like蛋白酶作用的肽類化合物結(jié)構(gòu)和活性數(shù)據(jù)均來自于文獻(xiàn)26。這些分子的抑制活性范圍IC50較寬,將其值轉(zhuǎn)化pIC50(-logIC50)值作為定量構(gòu)效關(guān)系分析中的因變量,使數(shù)值均勻地分布在4.25-6.50之間。整個(gè)數(shù)據(jù)集以近似3:1的比例,被分為兩個(gè)部分。其中,33個(gè)分子結(jié)構(gòu)和活性作為訓(xùn)練集,用于模型的構(gòu)建;而其余10個(gè)分子則加入到驗(yàn)證集中,用于測(cè)試該模型的準(zhǔn)確性。訓(xùn)練集和驗(yàn)證集的選取應(yīng)遵循隨機(jī)挑選的原則,在盡可能滿足更好的結(jié)構(gòu)差異性與多樣性的同時(shí),使訓(xùn)練集與驗(yàn)證集的分子比較均勻地分散在整個(gè)數(shù)據(jù)庫范圍中。我們將所有的分子結(jié)構(gòu)及其對(duì)應(yīng)的生理活性pIC50實(shí)驗(yàn)值列于表S1(Supporting Information)中。其中,表S1還包含基于配體和受體的CoMFA與CoMSIA模型得到的pIC50預(yù)測(cè)值,表S1粗體字顯示的部分是具有最佳預(yù)測(cè)性能的模型的pIC50預(yù)測(cè)值。
2.2構(gòu)象優(yōu)化及分子疊合
在建立3D-QSAR模型中,分子疊合是十分重要的環(huán)節(jié),其直接決定了模型的優(yōu)劣,而在分子疊合前需要進(jìn)行構(gòu)象優(yōu)化27。分子三維構(gòu)象優(yōu)化的步驟如下:
(1)采用Chemoffice軟件構(gòu)建化合物分子的二維和三維結(jié)構(gòu);
(2)導(dǎo)入到SYBYL中進(jìn)行能量?jī)?yōu)化,每個(gè)分子上加載Gasteiger-Hückel電荷28;
(3)利用Tripos力場(chǎng)29作為能量最小化和構(gòu)象搜索的共軛方法,且收斂條件設(shè)置為0.00419 kJ· mol-1;
(4)采用Powell能量梯度獲得穩(wěn)定構(gòu)象,且最大優(yōu)化次數(shù)設(shè)置為1000次,梯度設(shè)置為2.095 kJ· mol-1·nm-1。
在整個(gè)數(shù)據(jù)集中,我們選擇活性最大的1號(hào)分子進(jìn)行疊合,應(yīng)用SYBYL軟件,在1號(hào)分子上劃定所有分子的分子骨架,用紫色線顯示,如圖1a所示。在研究中,使用了三種不同的疊合方式:疊合I是基于配體的疊合,得到疊合圖為圖1b;疊合II是基于受體完成的,將所有分子都對(duì)接于受體中,得到對(duì)接結(jié)果最高打分的構(gòu)象進(jìn)行疊合,分子疊合圖為圖1c;疊合III也是一種基于受體的疊合方法,其中所有的分子構(gòu)象均來自于疊合II,隨后經(jīng)歷疊合I的疊合過程,并得到疊合圖為圖1d。
圖1 模板分子結(jié)構(gòu)及三種分子疊合圖Fig.1 Structure of template molecule and three kinds of molecular alignment
2.3CoMFA及CoMSIA研究
CoMFA及CoMSIA都是應(yīng)用最廣泛的合理藥物設(shè)計(jì)方法。CoMFA基于的假設(shè)是配體與受體的相關(guān)作用是非共價(jià)性的相互作用,利用分子力學(xué)力場(chǎng)(位阻場(chǎng)和靜電場(chǎng))可以較好地解釋它們的相關(guān)作用30。相比之下,CoMSIA是以CoMFA為基礎(chǔ),并對(duì)其進(jìn)行了補(bǔ)充和優(yōu)化。CoMSIA將CoMFA中的位阻場(chǎng)和靜電場(chǎng)進(jìn)一步細(xì)分為五個(gè)場(chǎng),即多引入了疏水、氫鍵供體和氫鍵受體場(chǎng),所以更細(xì)致地描述了配體與受體間的相互作用。此外,CoMSIA改變了探針原子與配體之間作用能的計(jì)算公式,非CoMFA傳統(tǒng)的庫侖(Coulomb)和范德華(Lennard-Jones)函數(shù),而采用了更加平滑的高斯(Gaussian)函數(shù)來計(jì)算分子場(chǎng)的數(shù)據(jù),從而有效避免了原子位置發(fā)生異常和配體周圍不同探針位點(diǎn)上勢(shì)能的差異31。而且在CoMSIA中,分子力場(chǎng)能量于不同的探針位點(diǎn)間衰減很快,致使計(jì)算過程自動(dòng)收斂,無需再設(shè)定能量截?cái)嘀怠?/p>
在本文的研究中,我們采用了CoMFA和CoMSIA兩種定量構(gòu)效關(guān)系方法,用來間接地反映了化合物分子與受體相互作用過程中兩者之間的非鍵相互作用特征,并預(yù)測(cè)良好的三維定量構(gòu)效關(guān)系模型。分子疊合后,分子空間取向基本一致,然后用一個(gè)粒子探針在分子周圍的空間中游走,計(jì)算粒子探針與分子之間的相互作用,并記錄下在空間不同坐標(biāo)中它們相互作用的能量值,進(jìn)而獲得分子力場(chǎng)數(shù)據(jù)。在CoMFA建模方法中,用sp3雜化的C+為探針,探針?biāo)鶐щ姾蔀?1.0e,半徑為0.152 nm,默認(rèn)步長(zhǎng)為0.2 nm,對(duì)整個(gè)三維網(wǎng)格進(jìn)行搜索,且位阻場(chǎng)和靜電場(chǎng)效應(yīng)能的閾值分別設(shè)定為默認(rèn)值209.29 kJ·mol-1,同時(shí)其它參數(shù)均為默認(rèn)值。對(duì)于CoMSIA建模中,用到的三維網(wǎng)格和CoMFA相同,以sp3雜化的C+為探針,電荷為+1.0e,半徑為0.100 nm,疏水性為+1.0,氫鍵供體和受體強(qiáng)度均為+1.0。此外,構(gòu)建這兩種的QSAR模型時(shí),衰減因子均設(shè)定為默認(rèn)值0.3。
2.4PLS分析與驗(yàn)證
本文采用偏最小二乘回歸分析(PLS)來建立CoMFA/CoMSIA描述符和抑制劑pIC50值之間的關(guān)系。在PLS的建模中,為了避免過度擬合問題,必須要對(duì)其中每個(gè)連續(xù)組分的重要性進(jìn)行嚴(yán)格測(cè)試,然后測(cè)試會(huì)在組分不重要時(shí)停止。交叉驗(yàn)證是一種能測(cè)試這種重要性的實(shí)用且可靠的方法32。我們運(yùn)用的交叉驗(yàn)證方法是抽一法(leave-oneout),它是依次從訓(xùn)練集中抽取一個(gè)樣本,余下的n-1個(gè)樣本作為訓(xùn)練集構(gòu)建模型,并用該模型預(yù)測(cè)抽取出的樣本,直到每個(gè)分子都被取出預(yù)測(cè)一次,最終可以得到最佳組分?jǐn)?shù)(OPN),與能夠衡量所建模型的內(nèi)部預(yù)測(cè)能力的重要統(tǒng)計(jì)指標(biāo),交叉驗(yàn)證相關(guān)系數(shù)(Q2)。而利用所有訓(xùn)練集的分子結(jié)合OPN,通過非交叉驗(yàn)證分析計(jì)算可以得出非交叉驗(yàn)證相關(guān)系數(shù)(R2ncv)、統(tǒng)計(jì)值(F)、訓(xùn)練集的估計(jì)標(biāo)準(zhǔn)誤差(SEE)和驗(yàn)證集的最小標(biāo)準(zhǔn)預(yù)測(cè)偏差(SEP)33。此外,根據(jù)預(yù)測(cè)驗(yàn)證集的活性,我們還可以得到能夠評(píng)估模型的外部預(yù)測(cè)能力的重要參數(shù),驗(yàn)證集分子的預(yù)測(cè)值相關(guān)系數(shù)(R2pre)。因此采用PLS能在避免過擬合現(xiàn)象和消除大多數(shù)可變因素的同時(shí),建立大量的QSAR方程34,35。最后,我們將CoMFA/CoMSIA結(jié)果以等勢(shì)圖的形式直觀地表示出來,并進(jìn)行下一步的分析。
2.5分子對(duì)接
為揭示肽類抑制劑小分子和蛋白大分子之間的結(jié)合機(jī)制,以及確定抑制劑的生物活性構(gòu)型,我們從Protein Data Bank數(shù)據(jù)庫36中選擇蛋白靶點(diǎn)的三維晶體結(jié)構(gòu)PDB ID:4YOI26,基于遺傳算法的GOLD 5.0版本軟件37,進(jìn)行分子對(duì)接研究。分子對(duì)接是利用配體小分子與受體蛋白大分子的完整結(jié)構(gòu)信息,計(jì)算把配體放在受體活性位點(diǎn)的位置,然后按照能量匹配、幾何匹配以及化學(xué)環(huán)境匹配的原則來評(píng)價(jià)抑制劑和蛋白的結(jié)合能力,并預(yù)測(cè)它們之間最佳的結(jié)合模式38。在實(shí)際過程中,不僅配體小分子具有柔性易產(chǎn)生多種構(gòu)象,而且受體大分子也是柔性且構(gòu)象可發(fā)生變化。在兩者分子對(duì)接過程中,它們的構(gòu)象會(huì)變化以找到一個(gè)能達(dá)成局部穩(wěn)定的契合點(diǎn)。因此,研究人員可以通過分子對(duì)接確定柔性的配體與受體蛋白的結(jié)合位置和空間取向,以此來研究藥物的作用機(jī)制并指導(dǎo)新藥物的設(shè)計(jì)。
分子對(duì)接之前,先對(duì)蛋白質(zhì)模型進(jìn)行修改,只留下單鏈,并在蛋白質(zhì)結(jié)構(gòu)中加入極性氫,且刪除原有的配體與雜原子。我們?cè)贕OLD中將受體大分子中每一個(gè)原子坐標(biāo)1 nm空間內(nèi)的氨基酸殘基所組成的空間范圍設(shè)為可能與配體小分子的結(jié)合位點(diǎn)。我們將43個(gè)分子分別對(duì)接到潛在的結(jié)合位點(diǎn)中,每次對(duì)接都會(huì)產(chǎn)生5個(gè)左右可能的構(gòu)象。所有構(gòu)象都會(huì)經(jīng)過GOLD評(píng)估,對(duì)其適當(dāng)性進(jìn)行打分。打分的分值計(jì)算考慮到了配體分子內(nèi)氫鍵和拉力的作用,以及配體與受體之間相互的氫鍵和范德華力的貢獻(xiàn)39。在分子對(duì)接過程中,打分函數(shù)采用GoldScore算法,其他參數(shù)采用軟件中的默認(rèn)值。
2.6分子動(dòng)力學(xué)
為了進(jìn)一步檢驗(yàn)分子對(duì)接研究中結(jié)果的穩(wěn)定性,以及構(gòu)建一個(gè)能夠表現(xiàn)分子所在真實(shí)環(huán)境的蛋白質(zhì)模型,即需要考慮蛋白質(zhì)大分子的動(dòng)態(tài)靈活性,我們應(yīng)用GROMACS軟件包中的GROMOS96力場(chǎng)40,41對(duì)分子對(duì)接得到的復(fù)合物開展了詳盡的分子動(dòng)力學(xué)模擬。分子動(dòng)力學(xué)模擬是研究生物大分子體系的常用計(jì)算方法。該法根據(jù)牛頓經(jīng)典力學(xué)原理,模擬分子體系的運(yùn)動(dòng)狀態(tài),計(jì)算在相空間中所有原子的運(yùn)動(dòng)軌跡,接著在從各種狀態(tài)的分子體系構(gòu)成的系綜中,選取樣本計(jì)算構(gòu)型積分,且以此為基礎(chǔ),進(jìn)行計(jì)算該體系的熱力學(xué)量等其他宏觀性質(zhì)。
進(jìn)行動(dòng)力學(xué)模擬之前,需要對(duì)模擬體系進(jìn)行一系列預(yù)處理。我們?cè)赑RODRG 2.542中獲得配體-蛋白質(zhì)復(fù)合物分子的拓?fù)湮募?,隨后將初始復(fù)合物放置在一個(gè)邊長(zhǎng)為10 nm的周期性立方體晶格體中,其中保持晶格上每一點(diǎn)與蛋白質(zhì)的距離至少是1 nm41。為了保證體系的電中性,蛋白質(zhì)復(fù)合物和水中所有原子均被隨機(jī)地放置于晶格盒子中,并用簡(jiǎn)單點(diǎn)電荷水來填充模擬體系的其余部分43。詳細(xì)的動(dòng)力學(xué)模擬步驟如下:首先為了避免分子間高能碰撞的可能性,對(duì)整個(gè)系統(tǒng)進(jìn)行2500步最快下降法操作,令其能量逐步最小化,再進(jìn)行2500步共軛梯度優(yōu)化。接著,當(dāng)整個(gè)系統(tǒng)的溫度和壓力上升到300 K和101.325 kPa且處于每納米8.38 kJ·mol-1穩(wěn)定力時(shí),對(duì)常溫常壓NPT系綜進(jìn)行50 ps恒壓力與500 ps密度平衡操作。最后,在周期性邊界條件下,設(shè)定積分步長(zhǎng)為0.002 ps,在保持恒溫的條件下進(jìn)行50 ns動(dòng)力學(xué)模擬。在整個(gè)模擬過程中,我們分別采用了SHAKE算法44與PME方法45來固定含碳?xì)滏I和處理庫侖作用。
3.1 CoMFA與CoMSIA模擬統(tǒng)計(jì)量分析
在之前的研究工作中,我們分別運(yùn)用了基于配體和基于受體的三種疊合方式對(duì)數(shù)據(jù)集中43個(gè)肽類分子進(jìn)行疊合,生成了六個(gè)不同的疊合模型,并對(duì)CoMFA兩種力場(chǎng)(位阻場(chǎng)與靜電場(chǎng))和CoMSIA五種力場(chǎng)(位阻場(chǎng)、靜電場(chǎng)、疏水場(chǎng)、氫鍵供體場(chǎng)和氫鍵受體場(chǎng))的所有3種與31種排列組合進(jìn)行了計(jì)算。為了提高模型的可對(duì)比性,所有的CoMFA和CoMSIA模型在建模過程中都使用相同的訓(xùn)練集(33個(gè)分子)和驗(yàn)證集(10個(gè)分子)。通過對(duì)比模型結(jié)果,得出基于配體得到的模型最優(yōu),因此我們主要分析基于配體的模型結(jié)果。這些模型的統(tǒng)計(jì)參數(shù)列于表1中。
表13D-QSAR結(jié)果匯總Table 1 Summary of 3D-QSAR results
對(duì)比得到了一個(gè)具有最佳預(yù)測(cè)性能的模型結(jié)果(表1中黑體字顯示),并用PLS法和抽一法對(duì)其進(jìn)行評(píng)估,獲得包括用來檢驗(yàn)所建模型是否有效的3個(gè)參數(shù)(Q2、OPN與R2ncv)和用來檢驗(yàn)?zāi)P皖A(yù)測(cè)能力的驗(yàn)證集R2pre。此外,還包括一些其他的統(tǒng)計(jì)參數(shù),如模型標(biāo)準(zhǔn)誤差,即F值、訓(xùn)練集的SEE和驗(yàn)證集的SEP,以及各個(gè)力場(chǎng)對(duì)模型的貢獻(xiàn)率。
從統(tǒng)計(jì)意義上來說,一個(gè)良好的3D-QSAR模型,前提必須滿足四個(gè)重要條件:Q2>0.5且,并有較低的SEE、較大的F值,才能很好地預(yù)示CoMFA和CoMSIA模型46。因此我們通過建模分析,最終以組合了擁有最高的的位阻場(chǎng)、靜電場(chǎng)、氫鍵供體場(chǎng)和氫鍵受體場(chǎng)的CoMSIA模型為最優(yōu)模型。其Q2=0.522、OPN= 9、R2ncv=0.996、F=602.246和SEE=0.042,可以看出這個(gè)模型的誤差較小且可靠性較高。當(dāng)被驗(yàn)證集驗(yàn)證時(shí),所得模型R2pre=0.904與SEP=0.440,展示出其良好的預(yù)測(cè)能力。所有這些數(shù)據(jù)結(jié)果都說明了我們所建立的基于配體疊合的CoMSIA模型是一個(gè)可靠的具有很好預(yù)測(cè)能力的模型。
同時(shí),圖2中顯示了在最優(yōu)的CoMSIA模型中,所有訓(xùn)練集分子和驗(yàn)證集分子的預(yù)測(cè)pIC50值和實(shí)際pIC50值的散點(diǎn)圖。其中訓(xùn)練集分子以方塊代表,而驗(yàn)證集分子以三角表示。從圖中可以看出,整個(gè)數(shù)據(jù)集都是密集地分布在回歸線的兩側(cè),而且抑制劑于模型中預(yù)測(cè)的pIC50值與實(shí)驗(yàn)中測(cè)得的pIC50值的差值較小。這表明訓(xùn)練集和驗(yàn)證集的密集度和誤差都很相似,未出現(xiàn)過擬合時(shí)訓(xùn)練集分布均勻,而驗(yàn)證集誤差較大(離回歸線很遠(yuǎn))的問題,且表明預(yù)測(cè)值和實(shí)際值的相關(guān)性良好,得到了比較可靠的構(gòu)效關(guān)系模型47,48。
圖2 基于配體的CoMSIA模型的預(yù)測(cè)與實(shí)際pIC50值的關(guān)聯(lián)圖Fig.2 Ligand-based correlation plot of the predicted versus the actual pIC50values based on the CoMSIAmodel
3.23D-QSAR等勢(shì)圖結(jié)果與分析
為進(jìn)一步形象化和具體化得到的CoMSIA模型中包含的信息和內(nèi)容,我們分析CoMSIA模型中四個(gè)力場(chǎng)的等勢(shì)圖。這些等勢(shì)圖能夠直觀地反映抑制劑小分子的生物活性與其化學(xué)結(jié)構(gòu)上的相互聯(lián)系,從而探究出對(duì)藥物分子活性有關(guān)的重要因素,進(jìn)而在實(shí)際藥物設(shè)計(jì)中提供指導(dǎo)49。在繪制等勢(shì)圖過程中,為了能夠更好地分析各種等勢(shì)圖的結(jié)果,一致將活性最高的分子,1號(hào)分子(pIC50= 6.481)選為模板分子,疊放于四張等勢(shì)圖(圖3)中,以便于對(duì)比。在該1號(hào)分子的三維結(jié)構(gòu)(圖4)中,R1、R2、R3和R4取代基是數(shù)據(jù)集內(nèi)其他分子結(jié)構(gòu)與1號(hào)分子的差別所在。此外,在該CoMSIA模型四個(gè)力場(chǎng)的等勢(shì)圖中,積極影響區(qū)域和消極影響區(qū)域貢獻(xiàn)比率分別設(shè)置為默認(rèn)值80%和20%。
圖3 結(jié)合1號(hào)分子的CoMSIA模型等勢(shì)圖Fig.3 CoMSIAcontour maps in combination with compound 1
圖4 1號(hào)分子結(jié)構(gòu)Fig.4 Structure of compound 1
圖3(a)為基于1號(hào)分子CoMSIA位阻場(chǎng)的等勢(shì)圖。其中,綠色和黃色云團(tuán)分別表示位阻基團(tuán)對(duì)分子活性的有利和有害區(qū)域。由圖可看出,R3取代基遠(yuǎn)離分子骨架的部分和遠(yuǎn)處R1與R2取代基之間有很大的綠色區(qū),說明R3取代基遠(yuǎn)位以及R1與R2取代基相近部分存在的較大基團(tuán)有益于提高分子活性,這一點(diǎn)從數(shù)據(jù)集中得到了充分的驗(yàn)證。例如,24號(hào)分子(pIC50=5.658)與34號(hào)分子(pIC50= 4.954)對(duì)比,以及1號(hào)分子(pIC50=6.481)與13號(hào)分子(pIC50=5.509)對(duì)比,它們的差別只在對(duì)應(yīng)位置存在較大基團(tuán),使位阻變大,有益于提高抑制活性。此外,R2取代基靠近分子骨架苯環(huán)的部分和R4取代基周圍存在兩塊黃色云團(tuán),顯示了位阻較大基團(tuán)在這片區(qū)域?qū)Ψ肿踊钚云鹨种谱饔谩@?號(hào)分子(pIC50=6.387)和28號(hào)分子(pIC50=5.377)的R4取代基中,H原子的位阻明顯小于叔丁基的位阻,使得28號(hào)分子較2號(hào)分子活性低很多。
圖3(b)為基于1號(hào)分子CoMSIA靜電場(chǎng)的等勢(shì)圖。其中,藍(lán)色云團(tuán)代表正電性積極影響區(qū)域;而紅色云團(tuán)表示負(fù)電性積極影響區(qū)域。圖中R4取代基臨近分子骨架的部分有大片藍(lán)色區(qū)域,說明此處正電性取代基的存在對(duì)提高抑制活性有利,如在此處帶有正電性的叔丁基基團(tuán)的分子活性較大。而圖中R2取代基臨近分子骨架的部分有大片紅色區(qū)域,另外還有小片紅色區(qū)域在R3取代基靠近分子骨架氧原子的部分,說明負(fù)電性取代基在此處有利。圖可見,R2取代基臨近分子骨架部分和R1與R3取代基遠(yuǎn)離分子骨架部分附近有較大的青色區(qū)域,說明該區(qū)域有氫鍵供體作用的取代基對(duì)提高分子活性有利。例如,最高活性的1號(hào)分子在R2取代基處有N―取代基(其上擁有氫原子作為氫鍵供體),其抑制活性遠(yuǎn)遠(yuǎn)高于其它無氫鍵供體基團(tuán)或N―取代基(其上未擁有氫原子)的類似物分子,即是一很好的佐證。總的來說此處氫鍵供體基團(tuán)能提高抑制劑活性。同時(shí),R3與R4取代基的遠(yuǎn)位有兩個(gè)極小紫色云團(tuán),其對(duì)分子活性具有一定消極影響,但從大趨勢(shì)來看,基本可以忽略。
圖3(d)為基于1號(hào)分子CoMSIA氫鍵受體場(chǎng)的等勢(shì)圖。其中,紫色云團(tuán)代表更依賴氫鍵受體區(qū)域,氫鍵受體在此處更益于生物活性,而紅色云團(tuán)則相反。從圖可見,紫色與紅色云團(tuán)基本都同時(shí)出現(xiàn)在R2與R3取代基靠近分子骨架部分,但它們都處于非同側(cè),此外還有一塊紅色云團(tuán)在R2取代基離骨架的遠(yuǎn)位。這使分子具有一定的敏感性,說明R2與R3取代基的空間旋轉(zhuǎn)角度不同,導(dǎo)致其受氫鍵受體場(chǎng)的影響有巨大差異。當(dāng)R2與R3取代基能作為氫鍵的接收基團(tuán)時(shí),且接收基團(tuán)都恰好旋轉(zhuǎn)到空間中氫鍵受體有益活性區(qū)域時(shí),積極作用達(dá)到最大;反之,消極作用最大。例如,在基于9號(hào)分子(pIC50=5.699)的CoMSIA氫鍵受體場(chǎng)的等勢(shì)圖(圖5)中,可以看到,9號(hào)分子的R2取代基上具有一個(gè)S原子和兩個(gè)N原子可作為氫鍵受體,R2取代基的角度已基本達(dá)到使它們靠近紫色積極云團(tuán)而遠(yuǎn)離紅色云團(tuán);同樣,該現(xiàn)象也可從這個(gè)分子的R3取代基的三個(gè)N原子上看到。氫鍵受體場(chǎng)對(duì)9號(hào)分子的積極作用最大,這可能是證明該分子在數(shù)據(jù)集中具有較高活性的原因之一。
圖5 結(jié)合9號(hào)分子的CoMSIA模型氫鍵受體輪廓圖Fig.5 CoMSIAH-bond acceptor contour map in combination with compound 9
綜合以上3D-QSAR的研究結(jié)果,我們分析了影響分子抑制活性的各基團(tuán)作用域(圖6)。我們能夠得出概括性的關(guān)于肽類分子的結(jié)構(gòu)性優(yōu)化方面的結(jié)論:
(1)對(duì)于R1取代基,遠(yuǎn)離分子骨架的部分(靠近R2取代基的周圍),如圖6紅色區(qū)域內(nèi),大位阻基團(tuán)和氫鍵供體基團(tuán)有利于分子活性的提高;
(2)對(duì)于R2取代基,臨近分子骨架的部分,如圖6深藍(lán)色區(qū)域內(nèi),負(fù)電性基團(tuán)、氫鍵供體基團(tuán)和氫鍵受體基團(tuán)有利于分子活性的提高;此外,靠近R1取代基的周圍,如圖6淺藍(lán)色區(qū)域內(nèi),大位阻基團(tuán)有利于分子活性的提高;
(3)對(duì)于R3取代基,臨近分子骨架的部分,如圖6深綠色區(qū)域內(nèi),負(fù)電性基團(tuán)和氫鍵受體基團(tuán)有利于分子活性的提高;此外,遠(yuǎn)離分子骨架的部分,如圖6淺綠色區(qū)域內(nèi),大位阻基團(tuán)和氫鍵供體基團(tuán)有利于分子活性的提高;
(4)對(duì)于R4取代基,臨近分子骨架的部分,如圖6黃色區(qū)域內(nèi),正電性取代基有利于分子活性的提高。
圖6 分子抑制作用的交互特性Fig.6 Interaction features of compound impacting the inhibitory effect color online
3.3分子對(duì)接研究
分子對(duì)接是一種用于驗(yàn)證3D-QSAR模型準(zhǔn)確性和穩(wěn)定性的有效工具。它能幫助我們更好地預(yù)測(cè)候選藥物小分子與靶點(diǎn)蛋白之間可能的結(jié)合構(gòu)象,進(jìn)而預(yù)測(cè)小分子的結(jié)合力和生物活性。在本文的研究中,為了闡明所選肽類3C-like蛋白酶抑制劑分子能否調(diào)節(jié)靶點(diǎn)蛋白,同時(shí)為了進(jìn)一步形象化地展現(xiàn)它們之間的作用機(jī)制及結(jié)合構(gòu)象,我們對(duì)整個(gè)數(shù)據(jù)集中所有43個(gè)肽類分子都進(jìn)行了對(duì)接研究。
為了之后能與John等26已做的X射線晶體實(shí)驗(yàn)結(jié)果作對(duì)比分析,我們采用相同的條件(存在結(jié)晶水)的受體和配體分子做了分子對(duì)接。另外,我們也對(duì)比了不存在結(jié)晶水的受體和配體分子的對(duì)接結(jié)果。在它們的對(duì)接結(jié)果中,獲得得分最高的合理性打分,其平均分分別為81.519分與78.922分。從打分結(jié)果來看,與存在結(jié)晶水的對(duì)接結(jié)果相比,沒有結(jié)晶水的對(duì)接打分變低了2.597分,受體存在結(jié)晶水的對(duì)接效果確實(shí)比較好,驗(yàn)證了John等26實(shí)驗(yàn)條件的合理性。通過對(duì)存在結(jié)晶水的對(duì)接結(jié)果分析得到了結(jié)晶水HOH523與配體分子產(chǎn)生的氫鍵距離圖(圖7),其中它們的氫鍵距離范圍從0.234至0.334 nm。一般來說,形成氫鍵的幾何標(biāo)準(zhǔn)是距離小于0.35 nm,且它們之間所連成的角度必須大于120°,此外0.32-0.40 nm之間的距離可以稱為是弱氫鍵50。因此,HOH523能與大部分的配體(30/43)產(chǎn)生較強(qiáng)的氫鍵作用,表明了由于它們的氫鍵存在使配體更加穩(wěn)定的結(jié)合于空腔中,進(jìn)而影響到整體打分。最終,我們選擇1號(hào)配體分子與存在結(jié)晶水的受體大分子進(jìn)行分子對(duì)接獲得的分?jǐn)?shù)最高的、最合理的結(jié)果來進(jìn)行接下來的分析。
圖8為1號(hào)分子和原來蛋白質(zhì)晶體中結(jié)合的分子互相疊合的結(jié)果,它們分別呈橙色和綠色。從圖中可明顯看出,兩個(gè)分子在受體中的位置極其相似,這有力地證明了分子對(duì)接結(jié)果的準(zhǔn)確性,能夠再現(xiàn)3C-like蛋白酶與其抑制劑分子間的實(shí)驗(yàn)結(jié)合空腔。
在配體藥物分子與受體蛋白相結(jié)合時(shí),其主要作用一般是與其周圍的氨基酸形成的疏水、氫鍵等相關(guān)作用。圖9展示了最高活性的1號(hào)分子在其周圍0.45 nm內(nèi)的蛋白質(zhì)空腔中的氨基酸殘基的對(duì)接作用模式,包含的氨基酸從圖中可以看出,抑制劑分子穩(wěn)定在蛋白質(zhì)空腔內(nèi),形成一個(gè)穩(wěn)定的空間區(qū)域。從圖中配體分子周圍的氨基酸類型判斷,該分子結(jié)合的空腔有兩個(gè)疏水位點(diǎn)(文中標(biāo)為S1和S2)。S1位點(diǎn)由氨基酸Ser24、Met25、ILe42、Pro45、Ala46以及Leu49等氨基酸構(gòu)成。1號(hào)抑制劑分子的R2取代基部分正結(jié)合于此。S1位點(diǎn)中的氨基酸大部分都具有脂肪族氨基酸的特性,它們共同促進(jìn)該疏水位點(diǎn)的形成。同時(shí),1號(hào)抑制劑分子的R3取代基部分也被一個(gè)疏水性的S2位點(diǎn)所固定。這些氨基酸包括Leu27、Phe143、Leu144、Cys145、Met165、His166、Met168以及Glu169等,它們也有較強(qiáng)的脂肪特性使得配體分子能夠緊密固定在這個(gè)空腔中。
圖7 結(jié)晶水HOH523與配體分子的氫鍵距離Fig.7 H-bond lengths between HOH523 and ligand molecules
此外,以上對(duì)接結(jié)果圖顯示,配體還通過不同方向的四個(gè)氫鍵緊密地結(jié)合在受體空腔內(nèi)。四個(gè)氫鍵分別為:
(1)配體分子R2取代基上的氮原子與結(jié)晶水HOH523氧原子形成的氫鍵(N―H…O,0.269 nm);
(2)配體分子R3取代基上的遠(yuǎn)離骨架部分的氮原子與氨基酸His166側(cè)鏈上氮原子形成的氫鍵(N…H―N,0.290 nm);
(3)配體分子骨架上的氧原子與氨基酸Glu169骨架上氮原子形成的弱氫鍵(O…H―N,0.359 nm);
(4)配體分子R3取代基上的靠近骨架部分的氮原子與氨基酸Glu169骨架上氮原子形成的弱氫鍵(N…H―N,0.352 nm)。
圖8 對(duì)接在3C-like蛋白酶中的1號(hào)分子的結(jié)合模式(橙色)和原1號(hào)分子(綠色)Fig.8 Binding mode of compound 1(orange)and the original molecule(green)docked in 3C-like protease
圖9 對(duì)接在3C-like蛋白酶中的抑制劑分子的結(jié)合模式Fig.9 Binding mode of inhibitor compound docked in 3C-like protease
與之前的CoMSIA等勢(shì)圖的結(jié)果結(jié)合分析,由CoMSIA氫鍵供體場(chǎng)的等勢(shì)圖(圖3(c))表明,在氫鍵供體場(chǎng)的積極作用域內(nèi),R2取代基上的氮原子作為氫鍵供體和結(jié)晶水HOH523形成氫鍵;而由CoMSIA氫鍵受體場(chǎng)的等勢(shì)圖(圖3(d))表明,在氫鍵受體場(chǎng)的積極作用域內(nèi),骨架上的氧原子和R3取代基上的靠近骨架部分的氮原子都作為氫鍵受體和氨基酸Glu169形成氫鍵。并且,CoMSIA靜電場(chǎng)的等勢(shì)圖(圖3(b))顯示,R4與R3取代基周圍分別存在有正電與負(fù)電云團(tuán),從對(duì)接結(jié)果中得知,在這兩個(gè)位置分別看到了正電性的氨基酸His41與Lys191和負(fù)電性的氨基酸Glu169。此外,由于Phe143、His166及His175等氨基酸中相對(duì)較大的苯環(huán)或雜環(huán)結(jié)構(gòu)的存在,也使得在R3取代基附近引入大型取代基會(huì)產(chǎn)生較強(qiáng)的位阻效應(yīng),這與CoMSIA位阻場(chǎng)的等勢(shì)圖(圖3(a))中大塊的綠色云圖所代表的含義也是十分一致的。這些區(qū)域中氫鍵的形成以及正負(fù)電性與含苯環(huán)或雜環(huán)結(jié)構(gòu)氨基酸的存在有利于提高生物的抑制活性。分子對(duì)接的結(jié)果與3D-QSAR等勢(shì)圖結(jié)果非常吻合,它們互為驗(yàn)證和補(bǔ)充,充分說明了3D-QSAR模型的合理性。
同時(shí),我們也將分子對(duì)接的結(jié)果與John等26已做的X射線晶體實(shí)驗(yàn)結(jié)果(PDB ID:4YOI、4YOG與4YOJ)進(jìn)行了對(duì)比。其中,氫鍵的對(duì)比結(jié)果如表2所示。從表中可以看出對(duì)接結(jié)果所得的四個(gè)氫鍵同樣在該實(shí)驗(yàn)結(jié)果中出現(xiàn),且它們的差別僅在作用力大小的不同,充分說明對(duì)接結(jié)果的正確性和穩(wěn)定性。此外,他們實(shí)驗(yàn)結(jié)果中的關(guān)鍵氨基酸,如氨基酸Met25、His41、Cys44、Ala46、Tyr54、Cys148、His166、Glu169和Gln192都出現(xiàn)在了分子對(duì)接的結(jié)合空腔中,也說明了分子對(duì)接是合理的,表明對(duì)接結(jié)果有很強(qiáng)的可靠性。
3.4分子動(dòng)力學(xué)研究
對(duì)于配體分子結(jié)合機(jī)制的探明是一個(gè)重要的研究步驟,因此我們需要獲得肽類抑制劑小分子和3C-like蛋白酶大分子復(fù)合物更加真實(shí)的結(jié)合機(jī)制。分子動(dòng)力學(xué)模擬正是能夠?qū)崿F(xiàn)這種優(yōu)化的方法之一51。分子動(dòng)力學(xué)模擬可以有效地表達(dá)分子體系的狀態(tài)和行為隨時(shí)間的變化情況。相比屬于半柔性模擬方法的分子對(duì)接研究,分子動(dòng)力學(xué)方法屬于動(dòng)態(tài)模擬的過程,能夠高效迅速地搜索出配體分子的低能量構(gòu)象,并對(duì)其構(gòu)象變化的軌跡進(jìn)行追蹤和表現(xiàn)。實(shí)驗(yàn)中,我們采用了水環(huán)境下分子動(dòng)力學(xué)模擬的方法來估計(jì)所得配體分子的結(jié)合親和性,并進(jìn)一步評(píng)估分子對(duì)接中結(jié)合機(jī)制模型的可靠性。
3.4.11號(hào)分子的分子動(dòng)力學(xué)研究
我們以分子對(duì)接中的配體-蛋白質(zhì)復(fù)合物作為初始結(jié)構(gòu),對(duì)1號(hào)分子進(jìn)行了50 ns的分子動(dòng)力學(xué)研究,其發(fā)生構(gòu)象改變的動(dòng)力學(xué)圖像如圖10所示。其中,我們采取了通過對(duì)最初構(gòu)象的監(jiān)測(cè)從而對(duì)結(jié)構(gòu)差異性進(jìn)行幾何學(xué)測(cè)量的方法,均方根偏差(RMSD),來進(jìn)一步地確保取樣方法的合理性與復(fù)合物大分子的動(dòng)力學(xué)穩(wěn)定性。圖10a分別展現(xiàn)了1號(hào)復(fù)合物(黑色)、蛋白(紅色)與1號(hào)配體(藍(lán)色)的分子動(dòng)力學(xué)模擬的RMSD軌跡。圖中可以看出復(fù)合物與蛋白的RMSD軌跡極為接近,復(fù)合物的RMSD軌跡大部分被蛋白的RMSD軌跡覆蓋。它們的RMSD軌跡測(cè)量值范圍從0.010至0.050 nm,在25 ns時(shí)達(dá)到0.040 nm,且在此后的模擬過程中,該波動(dòng)一直都保持在0.040 nm上下。而配體的RMSD軌跡,測(cè)量值范圍從0.005至0.030 nm,在5 ns后波動(dòng)都穩(wěn)定在0.020 nm左右。這些結(jié)果表明該分子動(dòng)力學(xué)軌道具有良好的平衡性,并且系統(tǒng)中的分子對(duì)接復(fù)合物保持相對(duì)穩(wěn)定。因此,我們采用1號(hào)復(fù)合物最后10 ns的平均結(jié)構(gòu)來研究,相比于只采用單純的蛋白質(zhì)晶體結(jié)構(gòu)而言,會(huì)有更好的準(zhǔn)確性和穩(wěn)定性。同時(shí),之前分子對(duì)接中得到的復(fù)合物也放在一起加以比較(圖10b),分子動(dòng)力學(xué)和分子對(duì)接中得到的配體分子分別設(shè)為橙色和綠色的棍狀(圖10c)。從圖10b中我們可以直觀地看到,1號(hào)分子在分子動(dòng)力學(xué)模擬與分子對(duì)接研究中占據(jù)的結(jié)合位點(diǎn)相同,且在構(gòu)型構(gòu)象上沒有明顯的差異,這印證了分子對(duì)接模型的合理性。但也有一點(diǎn)構(gòu)象上的改變,仍不能忽視:在結(jié)合空腔中,分子動(dòng)力學(xué)中配體的R1和R2取代基與分子對(duì)接相比,扭轉(zhuǎn)了一些角度。考慮到可能由于以上的構(gòu)象變化而帶來的作用力的改變,我們對(duì)于分子動(dòng)力學(xué)研究中得到復(fù)合物的結(jié)合機(jī)制也進(jìn)行了相應(yīng)分析。直觀的作用力構(gòu)成如圖11所示,結(jié)合作用力主要包括疏水作用和四個(gè)氫鍵作用力。
表2 分子對(duì)接和實(shí)驗(yàn)結(jié)果的氫鍵距離Table 2 H-bond lengths from docking and experiment
從圖11中可以得出,1號(hào)分子的R2取代基部分固定在一個(gè)疏水的口袋之中,這個(gè)口袋主要由以下氨基酸殘基組成:Gly23、Ser24、Met25、Thr26、Leu49以及Ser147。而R3取代基部分也被一個(gè)疏水性的脂肪區(qū)域所包裹,組成這個(gè)區(qū)域的氨基酸主要包括 Phe143、Leu144、Cys145、Gly146、His166、Met168、Glu169以及Leu170。顯然,這兩個(gè)結(jié)合區(qū)域與分子對(duì)接中分析出的疏水位點(diǎn)S1和S2是相一致的。另一方面,分子動(dòng)力學(xué)中配體分子結(jié)合的作用力與之前揭示的分子對(duì)接中的作用力是較為相似的,結(jié)晶水與氨基酸Glu169的重要作用力依然存在,說明了水分子介導(dǎo)的關(guān)鍵作用以及氨基酸Glu169的保守性,在一定程度上支持了分子對(duì)接過程中所得到的結(jié)論。但由于分子的R1和R2取代基扭轉(zhuǎn)了一些角度的原因使得分子與氨基酸His166的作用力消失,同樣與氨基酸Ser24產(chǎn)生了一個(gè)新的氫鍵作用力,使配體更有力地結(jié)合于對(duì)應(yīng)位點(diǎn)中,而這并不影響分子整體的結(jié)合位點(diǎn)。具體分子對(duì)接與分子動(dòng)力學(xué)氫鍵對(duì)比可由表3所示。總之,在水環(huán)境下的分子動(dòng)力學(xué)研究中,1號(hào)分子在活性位點(diǎn)處保持穩(wěn)定,所得的結(jié)果與分子對(duì)接結(jié)果也基本保持一致。
圖10 對(duì)接在3C-like蛋白酶中的1號(hào)分子的分子動(dòng)力學(xué)模擬結(jié)果Fig.10 MD-simulated results of compound 1 docked in 3C-like protease
圖11 在水環(huán)境下分子動(dòng)力學(xué)模擬的結(jié)合體的結(jié)合位點(diǎn)圖Fig.11 Plot of the in-water MD-simulated structures of the binding site
表3 分子對(duì)接與分子動(dòng)力學(xué)模擬的氫鍵分析Table 3 H-bond analysis from docking and MD simulation
3.4.2 2號(hào)分子的分子動(dòng)力學(xué)研究
除了1號(hào)分子的分子動(dòng)力學(xué)研究,我們還對(duì)2號(hào)分子進(jìn)行了50 ns分子動(dòng)力學(xué)研究。2號(hào)復(fù)合物(黑色)、蛋白(紅色)與2號(hào)配體(藍(lán)色)的RMSD圖以及結(jié)合模式圖如圖12所示。在圖12a所示的RMSD圖譜中,復(fù)合物與蛋白的RMSD數(shù)值在10 ns后穩(wěn)定在0.040 nm左右,而配體的數(shù)值則在20 ns后保持在0.020 nm上下,這顯示了一個(gè)穩(wěn)定的分子動(dòng)力學(xué)軌道。2號(hào)復(fù)合物最后10 ns的結(jié)合模式圖如圖12b所示,2號(hào)分子的R2取代基部分同樣固定在一個(gè)由Ser24、Met25、Thr26、Ala46、Asp47以及Leu49等氨基酸殘基組成的疏水口袋之中。而R3取代基部分也被一個(gè)由 Cys145、Gly146、Cys148、Met168、Glu169、Leu170以及His194等氨基酸組成的疏水脂肪區(qū)域所包裹。很顯然,這兩個(gè)小的疏水空腔與前文的位點(diǎn)S1和S2是基本相同的。這表明2號(hào)分子仍處于上述提及的疏水結(jié)合位點(diǎn)中。
圖12 對(duì)接在3C-like蛋白酶中的2號(hào)分子的分子動(dòng)力學(xué)模擬結(jié)果Fig.12 MD-simulated results of compound 2 docked in 3C-like protease
除了疏水作用,2號(hào)分子的分子動(dòng)力學(xué)所得到的氫鍵作用列于表4中,其中也列出了1號(hào)分子的分子動(dòng)力學(xué)中的氫鍵作用。對(duì)比表4的數(shù)據(jù)可以發(fā)現(xiàn),它們各自四個(gè)氫鍵作用中有三個(gè)氫鍵相同,并且可以發(fā)現(xiàn)2號(hào)分子與Gln192殘基又有一個(gè)新的氫鍵作用。分子在結(jié)構(gòu)方面的特征最終決定了它們各自動(dòng)力學(xué)模擬時(shí)的特殊構(gòu)造要求52。我們通過對(duì)比1號(hào)與2號(hào)分子結(jié)構(gòu)之間的差異進(jìn)行分析,它們的結(jié)構(gòu)特征圖如圖13所示。由圖13可以看出,它們的R1、R3與R4取代基以及R2取代基靠近分子骨架部分擁有共同的結(jié)構(gòu),而它們的差別僅于R2取代基遠(yuǎn)位部分,1號(hào)分子被噻吩環(huán)取代,而2號(hào)分子被苯環(huán)取代。大部分共同的結(jié)構(gòu)使它們處于位阻、靜電、氫鍵供體與氫鍵受體場(chǎng)的積極區(qū)域有諸多相同,但是由于噻吩環(huán)中具有一個(gè)S原子可作為氫鍵受體益于分子活性,導(dǎo)致1號(hào)分子活性略高于2號(hào)分子。因此,以上1號(hào)與2號(hào)分子結(jié)構(gòu)之間的微小差異導(dǎo)致了結(jié)合模式中氫鍵的細(xì)微改變,體現(xiàn)了R2取代基遠(yuǎn)位部分在配體與蛋白結(jié)合過程中的重要性。綜上所述,在1號(hào)與2號(hào)分子的分子動(dòng)力學(xué)研究中,配體分子在活性位點(diǎn)都保持穩(wěn)定,所得的兩個(gè)結(jié)果也基本保持一致。
表4 1號(hào)與2號(hào)分子的分子動(dòng)力學(xué)模擬的氫鍵分析Table 4 H-bond analysis from MD simulation of compound 1 and compound 2
圖13 1號(hào)與2號(hào)分子的結(jié)構(gòu)特征圖Fig.13 Structural features of compound 1 and compound 2
3.5新肽類化合物的設(shè)計(jì)與活性預(yù)測(cè)
基于43個(gè)肽類抑制劑對(duì)受體3C-like蛋白酶的抑制活性,建立了具有較強(qiáng)可靠性和良好預(yù)測(cè)能力的3D-QSAR模型。通過對(duì)模型的分析揭示了這類抑制劑分子結(jié)構(gòu)上的特點(diǎn),為更好地理解配體抑制劑和受體蛋白的作用機(jī)理提供了有效的幫助,同時(shí)為進(jìn)一步優(yōu)化此類抑制劑提供了可靠的依據(jù)。圖14顯示的是我們借助上述研究所獲得的基于3C-like蛋白酶上影響肽類分子抑制作用的關(guān)鍵性結(jié)構(gòu)特征。該圖將有助于指導(dǎo)篩選和開發(fā)更優(yōu)良的抗MERS-CoV藥物。
圖14 基于3C-like蛋白酶受體上肽類配體分子的交互特性Fig.14 Interaction features of peptidomimetic ligand molecule with the 3C-like protease receptor
圖15 新設(shè)計(jì)分子的結(jié)構(gòu)Fig.15 Structures of newly designed molecules
例如,以1號(hào)分子為模板,我們根據(jù)以上肽類分子的交互特性圖,通過修飾其各個(gè)取代基,設(shè)計(jì)了8個(gè)新的肽類3C-like蛋白酶抑制劑。這8個(gè)分子的分子結(jié)構(gòu)如圖15所示,CoMSIA模型預(yù)測(cè)它們的pIC50值均高于模板分子(6.48),關(guān)于其具體的藥物活性還需通過實(shí)驗(yàn)進(jìn)一步驗(yàn)證。
以計(jì)算機(jī)輔助藥物設(shè)計(jì)方法的理論和手段,借助3D-QSAR建模、分子對(duì)接和分子動(dòng)力學(xué)方法對(duì)一系列新合成的肽類3C-like蛋白酶抑制劑進(jìn)行了分子構(gòu)效關(guān)系以及配受體結(jié)合作用的研究。本文具體研究結(jié)論如下:
(1)建立了基于配體疊合的3D-QSAR模型,模型表現(xiàn)出較好的內(nèi)部一致性,并得到了CoMSIA模型等統(tǒng)計(jì)指標(biāo);
(2)對(duì)于R1取代基,遠(yuǎn)離分子骨架的部分(靠近R2取代基的周圍),大位阻基團(tuán)和氫鍵供體基團(tuán)有利于分子活性的提高;
(3)對(duì)于R2取代基,臨近分子骨架的部分,負(fù)電性基團(tuán)、氫鍵供體基團(tuán)和氫鍵受體基團(tuán)有利于分子活性的提高;同樣,靠近R1取代基的周圍,大位阻基團(tuán)也有利于分子活性的提高;
(4)對(duì)于R3取代基,臨近分子骨架的部分,負(fù)電性基團(tuán)和氫鍵受體基團(tuán)有利于分子活性的提高;同樣,遠(yuǎn)離分子骨架的部分,大位阻基團(tuán)和氫鍵供體基團(tuán)也有利于分子活性的提高;
(5)對(duì)于R4取代基,臨近分子骨架的部分,正電性取代基有利于分子活性的提高;
(6)分子對(duì)接結(jié)果顯示,抑制劑分子通過疏水作用以及與結(jié)晶水HOH523、氨基酸His166和Glu169產(chǎn)生的四個(gè)氫鍵作用使小分子穩(wěn)定在受體的空腔中,說明疏水性以及結(jié)晶水、氨基酸His166和Glu169在配體和受體結(jié)合過程中產(chǎn)生重要的作用;
(7)分子動(dòng)力學(xué)結(jié)果不僅與分子對(duì)接結(jié)果有較一致的結(jié)合空腔,證明了分子對(duì)接結(jié)果的可靠性,而且在結(jié)合空腔中,發(fā)現(xiàn)了兩個(gè)新的關(guān)鍵氨基酸Ser24與Gln192,它們與配體產(chǎn)生了兩個(gè)較強(qiáng)的氫鍵。
結(jié)論表明我們的3D-QSAR、分子對(duì)接和分子動(dòng)力學(xué)的研究得到了一個(gè)令人滿意的結(jié)果。此外,根據(jù)這些結(jié)論,一些新的具有潛在抑制活性的肽類化合物作為3C-like蛋白酶抑制劑被獲得。通過以上結(jié)論不但能夠深入了解肽類抑制劑的結(jié)構(gòu)特征以及與3C-like蛋白酶的結(jié)合作用,也為今后它的實(shí)驗(yàn)設(shè)計(jì)以及抗MERS-CoV藥物合成提供了新的方向。
Supporting Information:available free of charge via the internet at http://www.whxb.pku.edu.cn.
(1) Perlman,S.;Netland,J.Nat.Rev.Microbiol.2009,7,439. doi:10.1038/nrmicro2147
(2) Woo,P.C.;Huang,Y.;Lau,S.K.;Yuen,K.Y.Viruses 2010,2, 1804.doi:10.3390/v2081803
(3) Zaki,A.M.;Van Boheemen,S.;Bestebroer,T.M.;Osterhaus, A.D.;Fouchier,R.A.N.Engl.J.Med.2012,367,1814. doi:10.1056/NEJMoa1211721
(4) Rasmussen,S.A.;Gerber,S.I.;Swerdlow,D.L.Clin.Infect. Dis.2015,60,1686.doi:10.1093/cid/civ118
(5) Lau,S.K.;Li,K.S.;Tsang,A.K.;Lam,C.S.;Ahmed,S.; Chen,H.;Chan,K.H.;Woo,P.C.;Yuen,K.Y.J.Virol.2013, 87,8638.doi:10.1128/JVI.01055-13
(6) Chan,J.F.W.;Lau,S.K.P.;Woo,P.C.Y.J.Formos.Med. Assoc.2013,112,372.doi:10.1016/j.jfma.2013.05.010
(7) de Groot,R.J.;Baker,S.C.;Baric,R.S.;Brown,C.S.; Drosten,C.;Enjuanes,L.;Fouchier,R.A.;Galiano,M.; Gorbalenya,A.E.;Memish,Z.A.J.Virol.2013,87,7790. doi:10.1128/JVI.01244-13
(8) Haagmans,B.L.;Al Dhahiry,S.H.;Reusken,C.B.;Raj,V.S.; Galiano,M.;Myers,R.;Godeke,G.J.;Jonges,M.;Farag,E.; Diab,A.Lancet.Infect.Dis.2014,14,140.doi:10.1016/S1473-3099(13)70690-X
(9) Niemeyer,D.;Zillinger,T.;Muth,D.;Zielecki,F.;Horvath,G.; Suliman,T.;Barchet,W.;Weber,F.;Drosten,C.;Müller,M.A. J.Virol.2013,87,12489.doi:10.1128/JVI.01845-13
(10) Matthews,K.L.;Coleman,C.M.;van der Meer,Y.;Snijder,E. J.;Frieman,M.B.J.Gen.Virol.2014,95,874.doi:10.1099/ vir.0.062059-0
(11) Wang,Q.;Qi,J.;Yuan,Y.;Xuan,Y.;Han,P.;Wan,Y.;Ji,W.;Li, Y.;Wu,Y.;Wang,J.Cell Host Microbe 2014,16,328. doi:10.1016/j.chom.2014.08.009
(12) Woo,P.C.;Lau,S.K.;Li,K.S.;Tsang,A.K.;Yuen,K.Y. Emerg.Microbes.Infect.2012,1,e35.doi:10.1038/emi.2012.45
(13) Ziebuhr,J.;Snijder,E.J.;Gorbalenya,A.E.J.Gen.Virol.2000, 81,853.doi:10.1099/0022-1317-81-4-853
(14) Thiel,V.;Ivanov,K.A.;Putics,A.;Hertzig,T.;Schelle,B.; Bayer,S.;Wei?brich,B.;Snijder,E.J.;Rabenau,H.;Doerr,H. W.J.Gen.Virol.2003,84,2305.doi:10.1099/vir.0.19424-0
(15) Liu,W.;Zhu,H.M.;Niu,G.J.;Shi,E.Z.;Chen,J.;Sun,B.; Chen,W.Q.;Zhou,H.G.;Yang,C.Bioorg.Med.Chem.2014, 22,292.doi:10.1016/j.bmc.2013.11.028
(16) Kuo,C.J.;Liang,P.H.ChemBioEng Rev.2015,2,118. doi:10.1002/cben.201400031
(17) Chen,S.;Chen,L.;Tan,J.;Chen,J.;Du,L.;Sun,T.;Shen,J.; Chen,K.;Jiang,H.;Shen,X.J.Biol.Chem.2005,280,164. doi:10.1074/jbc.M408211200
(18) Ramajayam,R.;Tan,K.P.;Liu,H.G.;Liang,P.H.Bioorg. Med.Chem.Lett.2010,20,3569.doi:10.1016/j. bmcl.2010.04.118
(19) Thanigaimalai,P.;Konno,S.;Yamamoto,T.;Koiwai,Y.; Taguchi,A.;Takayama,K.;Yakushiji,F.;Akaji,K.;Chen,S.E.; Naser-Tavakolian,A.Eur.J.Med.Chem.2013,68,372. doi:10.1016/j.ejmech.2013.07.037
(20) Kang,C.M.;Zhao,X.H.;Wang,X.Y.;Cheng,J.G.;Lü,Y.T. Acta Phys.-Chim.Sin.2013,29,431.[康從民,趙緒浩,王新宇,程家高,呂英濤.物理化學(xué)學(xué)報(bào),2013,29,431.]doi:10.3866/ PKU.WHXB201211151
(21) Zhang,S.Z.;Zheng,C.;Zhu,C.J.Acta Phys.-Chim.Sin.2015, 31,2395.[張淑貞,鄭 超,朱長(zhǎng)進(jìn).物理化學(xué)學(xué)報(bào),2015,31, 2395.]doi:10.3866/PKU.WHXB201510142
(22) Pillaiyar,T.;Manickam,M.;Jung,S.H.Med.Chem.2015,5, 361.doi:10.4172/2161-0444.1000287
(23) Ren,Z.;Yan,L.;Zhang,N.;Guo,Y.;Yang,C.;Lou,Z.;Rao,Z. Protein Cell 2013,4,248.doi:10.1007/s13238-013-2841-3
(24) Deng,X.;StJohn,S.E.;Osswald,H.L.;O'Brien,A.;Banach,B. S.;Sleeman,K.;Ghosh,A.K.;Mesecar,A.D.;Baker,S.C.J. Virol.2014,88,11886.doi:10.1128/JVI.01528-14
(25) Tomar,S.;Johnston,M.L.;John,S.E.S.;Osswald,H.L.; Nyalapatla,P.R.;Paul,L.N.;Ghosh,A.K.;Denison,M.R.; Mesecar,A.D.J.Biol.Chem.2015,290,19403.doi:10.1074/ jbc.M115.651463
(26) John,S.E.S.;Tomar,S.;Stauffer,S.R.;Mesecar,A.D.Bioorg. Med.Chem.2015,23,6036.doi:10.1016/j.bmc.2015.06.039
(27) AbdulHameed,M.D.M.;Hamza,A.;Liu,J.J.;Zhan,C.G.J. Chem.Inf.Model.2008,48,1760.doi:10.1021/ci800147v
(28) Gasteiger,J.;Marsili,M.Tetrahedron 1980,36,3219. doi:10.1016/0040-4020(80)80168-2
(29) Clark,M.;Cramer,R.D.;Van Opdenbosch,N.J.Comput. Chem.1989,10,982.doi:10.1002/jcc.540100804
(30) Cramer,R.D.;Patterson,D.E.;Bunce,J.D.J.Am.Chem.Soc. 1988,110,5959.doi:10.1021/ja00226a005
(31) Klebe,G.;Abraham,U.;Mietzner,T.J.Med.Chem.1994,37, 4130.doi:10.1021/jm00050a010
(32) Edraki,N.;Das,U.;Hemateenejad,B.;Dimmock,J.R.Iran.J. Pharm.Res.2016,15,425.
(33) Li,X.L.;Ye,L.;Wang,X.X.;Wang,X.Z.;Liu,H.L.;Qian,X. P.;Zhu,Y.L.;Yu,H.X.Sci.Total Environ.2012,441,230. doi:10.1016/j.scitotenv.2012.08.072
(34) Shah,P.;Saquib,M.;Sharma,S.;Husain,I.;Sharma,S.K.; Singh,V.;Srivastava,R.;Shaw,A.K.;Siddiqi,M.I.Bioorg. Chem.2015,59,91.doi:10.1016/j.bioorg.2015.02.001
(35) Zhang,S.;Hou,B.;Yang,H.;Zuo,Z.Arch.Pharm.Res.2016, 1.doi:10.1007/s12272-016-0709-9
(36) Berman,H.M.;Battistuz,T.;Bhat,T.N.;Bluhm,W.F.;Bourne, P.E.;Burkhardt,K.;Iype,L.;Jain,S.;Fagan,P.;Marvin,J.; Padilla,D.;Ravichandran,V.;Schneider,B.;Thanki,N.; Weissig,H.;Westbrook,J.D.;Zardecki,C.Acta Crystallogr.D 2002,58,899.doi:10.1107/S0907444902003451
(37) Verdonk,M.L.;Cole,J.C.;Hartshorn,M.J.;Murray,C.W.; Taylor,R.D.Proteins 2003,52,609.doi:10.1002/prot.10465
(38) Duan,A.X.;Chen,J.;Liu,H.D.;Liu,X.H.;Lu,X.Q.J.Mol. Sci.2009,25,473.[段愛霞,陳 晶,劉宏德,劉秀輝,盧小泉.分析科學(xué)學(xué)報(bào),2009,25,473.]
(39) Arooj,M.;Sakkiah,S.;Kim,S.;Arulalapperumal,V.;Lee,K. W.PloS One 2013,8,e63030.doi:10.1371/journal. pone.0063030
(40) Van der Spoel,D.;Lindahl,E.;Hess,B.;Groenhof,G.;Mark,A. E.;Berendsen,H.J.C.J.Comput.Chem.2005,26,1701. doi:10.1002/jcc.20291
(41) Lindahl,E.;Hess,B.;Van Der Spoel,D.J.Mol.Model.2001,7, 306.doi:10.1007/s008940100045
(42) vanAalten,D.M.F.;Bywater,R.;Findlay,J.B.C.;Hendlich, M.;Hooft,R.W.W.;Vriend,G.J.Comput.Aid.Mol.Des. 1996,10,255.doi:10.1007/BF00355047
(43) Berendsen,H.J.C.;Postma,J.P.M.;van Gunsteren,W.F.; Hermans,J.Interaction Models for Water in Relation to Protein Hydration.In Intermolecular Forces;Springer Netherlands: Berlin,1981;pp 331-342.
(44) Ryckaert,J.P.;Ciccotti,G.;Berendsen,H.J.C.J.Comput. Phys.1977,23,327.doi:10.1016/0021-9991(77)90098-5
(45) Essmann,U.;Perera,L.;Berkowitz,M.L.;Darden,T.;Lee,H.; Pedersen,L.G.J.Chem.Phys.1995,103,8577.doi:10.1063/ 1.470117
(46) Alexander,G.;Alexander,T.J.Mol.Graph.Model.2002,20, 269.doi:10.1016/S1093-3263(01)00123-1
(47) Kamsri,P.;Punkvang,A.;Hannongbua,S.;Saparpakorn,P.; Pungpo,P.RSC Adv.2015,5,52926.doi:10.1039/C5RA08103C
(48) Gao,X.;Han,L.;Ren,Y.Molecules 2016,21,591.doi:10.3390/ molecules21050591
(49) Da,C.X.;Mooberry,S.L.;Gupton,J.T.;Kellogg,G.E. J.Med.Chem.2013,56,7382.doi:10.1021/jm400954h
(50) Jeffrey,G.A.An Introduction to Hydrogen Bonding.In Topics in Physical Chemistry;Oxford University Press:NewYork,1997.
(51) Minini,L.;Alvarez,G.;González,M.;Cerecetto,H.;Merlino, A.J.Mol.Graph.Model.2015,58,40.doi:10.1016/j. jmgm.2015.02.002
(52) Wang,J.;Li,F.;Li,Y.;Yang,Y.;Zhang,S.;Yang,L.Mol. BioSyst.2013,9,2296.doi:10.1039/c3mb70105k
QSAR,Molecular Docking and Molecular Dynamics of 3C-like Protease Inhibitors
LIN Feng1FU Xin-Mei2,*WANG Chao1JIANG Si-Yu1WANG Jing-Hui1ZHANG Shu-Wei1YANG Ling3LIYan1,*
(1School of Chemical Engineering,Dalian University of Technology,Dalian 116024,Liaoning Province,P.R.China;2State Key Laboratory of Fine Chemicals,Dalian University of Technology,Dalian 116024,Liaoning Province,P.R.China;3Laboratory of Pharmaceutical Resource Discovery,Dalian Institute of Chemical Physics,Chinese Academy of Sciences, Dalian 116023,Liaoning Province,P.R.China)
3C-like protease is an extremely important protease involved in the multiplicative process of coronaviruses,includingthedeadlyMiddleEast respiratorysyndromecoronavirus(MERS-CoV).3C-likeprotease has become a hot research topic in the field of coronavirology.For the first time,a set of ligand-and receptorbased three-dimensional quantitative structure-activity relationships(3D-QSAR)models were carried out via comparative molecular field analysis(CoMFA)and comparative molecular similarity indices analysis(CoMSIA)to explore the structure-activity correlation of 43 peptidomimetic inhibitors of the 3C-like protease of the bat coronavirus HKU4(HKU4-CoV),which belongs to the same 2c lineage as MERS-CoV and shows high sequence similarity with MERS-CoV.Based on the ligand-based alignment,an optimal CoMSIAmodel(yielded by steric, electrostatic,H-bond donor and H-bond acceptor fields)was obtained with good predictive power of Q2=0.522, R2ncv=0.996 and R2pre=0.904(Q2:cross-validated correlation coefficient,R2ncv:non-cross-validated correlation coefficient,R2pre:predicted correlation coefficient for the test set of compounds).Molecular docking and molecular dynamics simulations were performed according to this model to further determine the interaction mechanism between ligands and the receptor.The experimental results show:(1)based on the optimal CoMSIAmodel,the 3D contour maps vividly illustrate that the molecular biological activity is influenced by the steric,electrostatic, H-bond donor and H-bond acceptor interactions of molecular groups.(2)Based on the docking analysis, hydrophobicity,crystal water,His166andGlu169haveimportantrolesintheligandsandreceptorbindingprocess. (3)Molecular dynamics(MD)simulations were carried out for further verification of the reliability of the docking model,and provide two new key residues,Ser24 and Gln192,which have two strong hydrogen bonds with the ligands.Some new compounds were obtained based on the modeling that are potential peptidomimetic inhibitors of 3C-like protease.These results help establish the binding mechanism between 3C-like protease and peptidomimetic inhibitors,and provide a valuable reference for future anti-MERS-CoV drug design.
MERS-CoV;3C-like protease;Peptidomimetic inhibitor;3D-QSAR;Molecular docking; Molecular dynamics
O641
10.3866/PKU.WHXB201608121
Received:May 31,2016;Revised:August 9,2016;Published online:August 12,2016.
*Corresponding authors.FU Xin-Mei,Email:fuxinmei@dlut.edu.cn;Tel:+86-411-84986206.LI Yan,Email:yanli@dlut.edu.cn; Tel:+86-411-84986062.
The project was supported by the National Natural Science Foundation of China(11201049).
國家自然科學(xué)基金(11201049)資助項(xiàng)目