孟現(xiàn)美　張少龍　張慶剛(山東師范大學(xué)物理與電子科學(xué)學(xué)院，濟(jì)南250014)

分子動(dòng)力學(xué)模擬別構(gòu)抑制劑Efavirenz對HIV-1逆轉(zhuǎn)錄酶的作用

孟現(xiàn)美張少龍*張慶剛*
(山東師范大學(xué)物理與電子科學(xué)學(xué)院，濟(jì)南250014)

為了理解非核苷類逆轉(zhuǎn)錄酶抑制劑(NNRTIs)與HIV-1逆轉(zhuǎn)錄酶(RT)的相互作用機(jī)制，利用新力場ff12SB對未結(jié)合和結(jié)合Efavirenz(EFV)逆轉(zhuǎn)錄酶的三種RT大分子體系分別進(jìn)行了100 ns的長時(shí)間動(dòng)力學(xué)模擬。通過分析EFV對RT結(jié)構(gòu)的影響、不同殘基柔性和不同體系構(gòu)象的動(dòng)力學(xué)行為等，發(fā)現(xiàn)EFV的結(jié)合會(huì)導(dǎo)致RT結(jié)構(gòu)變化，從而影響RT的活性；證實(shí)了EFV的“分子楔”作用；還發(fā)現(xiàn)EFV的結(jié)合不但引起“拇指關(guān)節(jié)炎”，而且引起輕度“手指關(guān)節(jié)炎”；整個(gè)模擬過程中沒有出現(xiàn)不同構(gòu)象間的躍遷，但是無別構(gòu)分子時(shí)的RT張開構(gòu)象表現(xiàn)出明顯的閉合傾向。這些結(jié)果有助于理解NNRTIs的抑制機(jī)制和RT構(gòu)象變化的動(dòng)力學(xué)性質(zhì)。另外，還比較分析了模擬方法對計(jì)算結(jié)果的影響，對大分子體系的動(dòng)力學(xué)模擬具有重要借鑒意義。

HIV-I逆轉(zhuǎn)錄酶；逆轉(zhuǎn)錄酶抑制劑；別構(gòu)抑制劑；分子動(dòng)力學(xué)模擬；構(gòu)象

doi:10.3866/PKU.WHXB201511302

1　引言

逆轉(zhuǎn)錄酶(RT)在I型人體免疫缺陷病毒(HIV-1)復(fù)制過程中，可以催化病毒RNA基因組使之逆轉(zhuǎn)錄為雙鏈DNA病毒前體，因此是發(fā)明抗艾滋病(anti-AIDS)新藥的一個(gè)重要靶酶1。目前，以RT為靶治療艾滋病的藥物中，非核苷類RT抑制劑(NNRTIs)被認(rèn)為是一類具有前景的特效藥2－4。NNRTIs也稱為別構(gòu)抑制劑，結(jié)合在RT的非活性位，可以避免活性位點(diǎn)或底物的選擇性壓力，具有廣譜抗HIV-1耐藥性的特點(diǎn)。如圖1A所示，RT是一個(gè)由p66和p51兩個(gè)亞基組成的異二聚體5。p66亞基由DNA聚合酶和RNase H結(jié)構(gòu)域構(gòu)成，其中，DNA聚合酶結(jié)構(gòu)域形狀類似于握杯的右手形，它的3個(gè)子域形象地稱之為手指(fingers)、手掌(palm)和拇指(thumb)子域，并通過所謂的連接子(connection)與RNase H結(jié)構(gòu)域連接在一起。p51亞基具有與p66聚合酶結(jié)構(gòu)域同樣的氨基酸序列，但折疊方式不同。如圖1B所示，在p66的拇指和手掌連接處有一個(gè)疏水性口袋結(jié)構(gòu)(NNRTI-BP)，NNRTIs可以進(jìn)入并束縛在袋內(nèi)6。圖中顯示的EFV (Efavirenz)屬于第二代NNRTIs類藥物。

圖1　逆轉(zhuǎn)錄酶(RT)和別構(gòu)抑制劑(EFV)的結(jié)構(gòu)Fig.1　Structures of reverse-transcriptase(RT)and the allosteric inhibitor(EFV)

然而，關(guān)于NNRTIs和RT的作用機(jī)制當(dāng)前還無一致的理解。一般觀點(diǎn)是建立在實(shí)驗(yàn)晶體數(shù)據(jù)分析基礎(chǔ)上的“拇指關(guān)節(jié)炎”或“分子楔”模型7－9，但最近的熒光實(shí)驗(yàn)卻意味NNRTIs的作用會(huì)使手指-拇指“握度松動(dòng)(grip loosening)”，從而抑制聚合催化過程10。理論上也相應(yīng)地開展了許多工作11－16，主要是利用分子動(dòng)力學(xué)方法研究具有不同初始結(jié)構(gòu)和不同底物的RT動(dòng)力學(xué)性質(zhì)。由于RT體系過大，所以模擬時(shí)間一般只有幾個(gè)納秒，所得計(jì)算結(jié)果也不同。例如，文獻(xiàn)11(3 ns)和12(2.5 ns)的模擬結(jié)果支持“拇指關(guān)節(jié)炎”模型；文獻(xiàn)13 (1 ns)模擬結(jié)果表明，去掉NNRTIs后拇指子域會(huì)迅速從張開塌向閉合狀態(tài)轉(zhuǎn)變；文獻(xiàn)14(1.1 ns)模擬結(jié)果表明，底物會(huì)增加RT的柔性。近來，兩個(gè)較長時(shí)間的動(dòng)力學(xué)模擬結(jié)果更是引人關(guān)注：Ivetac和McCammon15利用GROMACS程序和GROMOS力場，對具有不同原子初速度的體系副本進(jìn)行了30 ns的模擬，發(fā)現(xiàn)無論APO(無任何抑制劑)還是結(jié)合NNRTI的RT體系，其構(gòu)象發(fā)生轉(zhuǎn)變的幾率都相當(dāng)可觀(分別為50%和25%)；Wright等16用AMBER9和ff03力場進(jìn)行的計(jì)算(模擬時(shí)間達(dá)100 ns，但步長較大為4 fs)發(fā)現(xiàn)，被別構(gòu)抑制劑撐開的拇指子域可以重返類似APO的閉合狀態(tài)。顯然，這些模擬結(jié)果有悖于一般的“拇指關(guān)節(jié)炎”觀點(diǎn)，意味著RT的構(gòu)象變化和NNRTI的作用機(jī)制可能并不像一般認(rèn)為的那樣。換句話說，無論APO還是加上NNRTI后，RT體系的能量曲面上都可能存在兩個(gè)極值點(diǎn)，分別對應(yīng)張合狀態(tài)，且兩個(gè)極值點(diǎn)間過渡態(tài)的勢壘不高，常溫?zé)崃W(xué)運(yùn)動(dòng)就會(huì)引起從開放狀態(tài)到閉合狀態(tài)的轉(zhuǎn)變。如果情況果真如此，藥就將失去抑制功能，抗逆轉(zhuǎn)錄病毒治療必須考慮這種情況。因此，應(yīng)該弄清這些模擬結(jié)果的可靠性，或者說了解這些現(xiàn)象的發(fā)生幾率。對這些問題的闡明將有助于深入理解NNRTIs的抑制機(jī)制和聚合酶結(jié)構(gòu)域構(gòu)象變化的動(dòng)力學(xué)性質(zhì)，從而為設(shè)計(jì)更有效的抑制劑提供理論啟發(fā)。

分子動(dòng)力學(xué)模擬結(jié)果的可靠性依賴于模擬方法和策略，先進(jìn)的力場、較小的步長、較長的模擬時(shí)間以及體系動(dòng)力學(xué)過程的充分穩(wěn)定平衡將有助于改進(jìn)計(jì)算結(jié)果，得到更可信的信息17－19。新的GPU大規(guī)模并行計(jì)算技術(shù)的出現(xiàn)為生物大分子模擬計(jì)算提供了更加高效的平臺；新力場ff12SB20,21考慮到了側(cè)鏈χ1扭轉(zhuǎn)勢的重要性僅次于骨架扭轉(zhuǎn)勢，不只對蛋白的主鏈φ′/ψ′扭轉(zhuǎn)勢參數(shù)做了進(jìn)一步調(diào)整，并且也修正了17個(gè)氨基酸側(cè)鏈的扭轉(zhuǎn)勢參數(shù)，從而改善了勢函數(shù)的精度。

本文采用先進(jìn)的ff12SB力場，使用GPU計(jì)算技術(shù)，利用AMBER12程序包20，對結(jié)合和未結(jié)合別構(gòu)抑制劑EFV的三種RT體系進(jìn)行了模擬計(jì)算。為了敘述簡潔，結(jié)合EFV的RT體系簡稱EBR；未結(jié)合EFV，拇指子域?yàn)榇蜷_的RT體系簡稱為APOO，拇指子域?yàn)殚]合的簡稱為APOC。研究目的是，(1)別構(gòu)抑制劑EFV是如何結(jié)合到RT的非活性位并影響RT構(gòu)象；(2)EFV結(jié)合到非活性位后對RT穩(wěn)定性(即結(jié)構(gòu)柔性)會(huì)產(chǎn)生怎樣的影響；(3)通過先進(jìn)的力場、較小的時(shí)間步長、較長的模擬時(shí)間和嚴(yán)格的模擬初始結(jié)構(gòu)，重點(diǎn)探究不同初始構(gòu)象條件下，RT拇指子域構(gòu)象的動(dòng)力學(xué)演化性質(zhì)，以期驗(yàn)證張合兩個(gè)狀態(tài)間可能的轉(zhuǎn)變信息，進(jìn)而弄清NNRTIs和RT的作用機(jī)制。

2　理論計(jì)算方法

體系的初始構(gòu)象取自蛋白質(zhì)庫(PDB ID：1DLO，1IKW和1HQU)。1DLO中只包含RT，1IKW中包含RT和EFV，1HQU中包含RT和HBY (別構(gòu)抑制劑)。其中，所有帶別構(gòu)抑制劑的RT晶體結(jié)構(gòu)都是不完整的，在1IKW中，p51亞基缺失殘基217－231。由于1HQU的晶體結(jié)構(gòu)具有較高的分辨率(0.27 nm)，并結(jié)合了別構(gòu)抑制劑，所以我們借助該晶體補(bǔ)全1IKW缺失的殘基，補(bǔ)全的1IKW中，RT共擁有983個(gè)殘基。另外，相對于1DLO，1IKW中的RT殘基序列有兩處發(fā)生變異，即p66亞基上S280C和E478Q。對于1DLO和1IKW中的RT，因?yàn)榻M成RT的各氨基酸參數(shù)預(yù)存在庫文件中，所以直接利用AMBERTools12的LEaP模塊，根據(jù)ff12SB力場為RT賦鍵參數(shù)和靜電參數(shù)，并補(bǔ)全缺失的氫原子；對1IKW中的EFV，其力參數(shù)和靜電信息是沒有預(yù)存的，需要計(jì)算生成，我們利用sqm模塊中的AM1-BCC擬合EFV的電荷分布，其力場參數(shù)采用GAFF力場參數(shù)22和擬合電荷時(shí)在Antechamber程序中生成的參數(shù)。這樣，處理過的1IKW中的RT和EFV結(jié)合在一起組成體系EBR，RT單獨(dú)作為體系A(chǔ)POO；處理過的1DLO中的RT作為體系A(chǔ)POC。三個(gè)體系外圍都加了1 nm厚的截角八面體水盒子，并分別添加10、10和8個(gè)氯離子使之呈電中性。

圖2　三個(gè)NNRTI-BP體系中結(jié)合口袋入口處的初始構(gòu)象Fig.2　Starting conformations of entrance of the binding pocket in the three RT systems

為了消除原子間的過近接觸，采用AMBER12的PMEMD模塊對三個(gè)體系進(jìn)行能量最小化，先進(jìn)行25000步最陡下降法計(jì)算，后進(jìn)行25000步共軛梯度法計(jì)算，這樣就得到了具有不同初始構(gòu)象的三個(gè)模擬計(jì)算樣本。圖2示出了它們的別構(gòu)袋口處構(gòu)象，可以清楚看出，EBR、APOO和APOC的袋口分別處于張開、張開、閉合。接下來將上述三個(gè)體系的溫度在300 ps內(nèi)從0 K加熱到300 K，然后進(jìn)行200 ps的常溫常壓動(dòng)力學(xué)平衡計(jì)算，最后對三個(gè)體系分別進(jìn)行100 ns的分子動(dòng)力學(xué)模擬。計(jì)算中，時(shí)間步長取2 fs，每隔5 ps記錄一次坐標(biāo)文件，并使用SHAKE算法限制含氫原子的鍵的變化23。模擬過程中同時(shí)計(jì)算監(jiān)測系統(tǒng)的總能量、動(dòng)能和密度隨時(shí)間的變化情況。文中所有的結(jié)構(gòu)圖均采用VMD軟件24,25完成。

圖3EFV與RT結(jié)合后引起的RT結(jié)構(gòu)變化Fig.3　RT structure changes caused by binding EFV

3　結(jié)果與討論

3.1EFV與RT結(jié)合引起的結(jié)構(gòu)變化

NNRTI-BP由p66亞基上許多含芳環(huán)的殘基(Tyr181、Tyr188、Phe227、Trp229、Tyr232)，疏水殘基(Lys101、Lys103、Ser105、Val106、Val108、Asp192、Glu224等)以及p51亞基的Glu138組成?？诖娜肟谖挥趐66和p51亞基的接合部，由p66亞基的Leu100、Lys101、Lys103、Val179、Tyr181和p51亞基的Glu138環(huán)繞組成(圖2)7,26。通過圖2可以看出，沒有結(jié)合EFV的APOC體系，在NNRTIBP相應(yīng)位置存在一個(gè)表面凹陷，袋口是閉合的。當(dāng)EFV結(jié)合到RT上時(shí)，袋口張開，即RT結(jié)構(gòu)發(fā)生了變化，這種變化可以分為結(jié)構(gòu)上的短程和長程變化(圖3)。短程變化如圖3A所示，主要是構(gòu)成NNRTI-BP的氨基酸位置和取向的變化。變化比較明顯的是Tyr181(淺藍(lán)色)、Tyr188(大紅色)和Trp229(灰色)側(cè)鏈，未結(jié)合EFV時(shí)，這三個(gè)氨基酸的芳環(huán)側(cè)鏈指向疏水區(qū)核心，EFV的結(jié)合造成了側(cè)鏈顯著的扭轉(zhuǎn)，使之遠(yuǎn)離疏水核心，從而制造出一個(gè)可以容納小分子抑制劑的疏水空間。在這個(gè)過程中Val106(粉色)和Val108(紫紅色)兩個(gè)殘基發(fā)生了明顯的外向位移。進(jìn)入疏水口袋的EFV通過氫鍵與芳環(huán)之間的π鍵和van der Waals等作用與口袋周圍的氨基酸形成穩(wěn)定的復(fù)合物。為了研究EFV導(dǎo)致的RT長程結(jié)構(gòu)變化，將兩個(gè)體系的手指子域中Cα原子進(jìn)行疊合，如圖3B所示，即保持手指位置相對不變，然后測量其它子域方位的相對變化。發(fā)現(xiàn)，EFV的結(jié)合使拇指子域外旋約38°，連接子和RNase H結(jié)構(gòu)域的位置也分別遠(yuǎn)離手指12°和15°。另外我們還發(fā)現(xiàn)，RT結(jié)合EFV后，p51亞基上的手指、手掌、拇指子域和連接子也分別移動(dòng)9°、12°、14°和18°。

3.2RMSD變化

均方根偏差(RMSD)是衡量體系穩(wěn)定性的重要參量20,27，它反映蛋白質(zhì)分子在動(dòng)力學(xué)模擬過程中偏離初始結(jié)構(gòu)的程度。我們使用AMBER12中的ptraj程序分析了EBR、APOO和APOC三體系的骨架碳原子Cα相對于初始結(jié)構(gòu)的RMSD隨時(shí)間的變化(圖4)?？梢钥闯觯谡麄€(gè)模擬過程中，體系EBR(紅線)結(jié)構(gòu)的RMSD無論其數(shù)值大小還是漲落幅度均明顯低于APOO(黑線)和APOC(藍(lán)線)，說明體系EBR的穩(wěn)定性優(yōu)于APOO和APOC。

圖4　三體系骨架碳原子Cα的RMSD隨時(shí)間的變化Fig.4　Time dependences of root mean square deviation (RMSD)of the backbone Cαatoms for the three systems

具體看，EBR體系的RMSD變化情況比較簡單，最初0－12 ns不太穩(wěn)定，是一個(gè)由小變大趨于平衡的過程，12 ns之后直到模擬結(jié)束的100 ns，其值基本穩(wěn)定在0.24 nm附近，上下浮動(dòng)范圍約為0.04 nm。對于APOO和APOC兩個(gè)體系，在整個(gè)模擬時(shí)間內(nèi)，其RMSD表現(xiàn)出了相似的特征，在開始的0－40 ns內(nèi)，兩個(gè)體系的RMSD無論在變化趨勢還是波動(dòng)幅度上幾乎都是重合的，RMSD值由0.2 nm上升到0.4 nm；在之后的整個(gè)模擬時(shí)間內(nèi)，APOO和APOC的RMSD也表現(xiàn)了相似的波動(dòng)性，只是APOC的RMSD值較小波動(dòng)幅度較大，二者的RMSD平均值分別約為0.40和0.36 nm，浮動(dòng)范圍約為0.04和0.10 nm。

EBR體系較小的RMSD值表明，結(jié)合別構(gòu)抑制劑EFV后，RT結(jié)構(gòu)穩(wěn)定性變好，或者說EFV增加了體系的剛性，這正是一般認(rèn)為的別構(gòu)抑制劑NNRTIs的“楔子”作用，它阻止了相關(guān)子域的柔性。對于APOO和APOC兩個(gè)體系，其RMSD在整個(gè)模擬時(shí)間內(nèi)表現(xiàn)出的相似性意味著，沒有抑制劑時(shí)，張開狀態(tài)(APOO)與閉合狀態(tài)(APOC)具有類似的動(dòng)力學(xué)行為。值得注意的是，這時(shí)APOC的RMSD雖然具有較小的平均值但卻具有很寬的波動(dòng)幅度，其最大值與APOO的相近，而最小值卻與EBR的相當(dāng)，意味著APOC的狀態(tài)很不穩(wěn)定，結(jié)構(gòu)變化忽大忽小。而與之相對，EBR總是維持在結(jié)構(gòu)變化較小的狀態(tài)，APOO傾向于保持在結(jié)構(gòu)變化較大的狀態(tài)。這說明APOO和APOC兩體系比EBR具有較大結(jié)構(gòu)柔性，而APOC的結(jié)構(gòu)在動(dòng)力學(xué)過程中更易發(fā)生變化，即時(shí)柔時(shí)剛。顯然，不同體系所表現(xiàn)的這些不同動(dòng)力學(xué)行為只有通過較長的模擬時(shí)間才能獲得充分取樣，這解釋了為什么之前文獻(xiàn)13,14通過短時(shí)間模擬得到的結(jié)果是矛盾的。

3.3B因子分析

B因子(B-factor)描述體系中各個(gè)殘基在動(dòng)力學(xué)過程中偏離其平均結(jié)構(gòu)的程度，反映殘基的柔軟程度，可以用來研究分子動(dòng)力學(xué)過程中蛋白的一些性質(zhì)變化20。圖5給出了EBR、APOO和APOC三體系中所有殘基的B因子，它們?nèi)∽阅M過程中體系的RMSD數(shù)值較小且穩(wěn)定的時(shí)段，分別為，EBR：74－88 ns，APOO：12－32 ns，APOC：48－62 ns。

圖5A是三體系中RT全部983個(gè)殘基的B因子變化情況?？梢钥闯?，p51和p66亞基中距別構(gòu)抑制劑EFV結(jié)合位較遠(yuǎn)的連接子和RNase H結(jié)構(gòu)域，即320號之后的殘基，三體系相應(yīng)殘基的B因子變化不大，說明別構(gòu)抑制劑的加入對這些部位影響不大。然而，距離別構(gòu)位較近的聚合酶活性區(qū)，包括手指、手掌和拇指子域，即1到320號殘基，三體系相應(yīng)殘基的B因子差別非常明顯。

拇指和手指子域的構(gòu)象變化和柔性直接關(guān)系到RT的催化功能，所以有必要詳細(xì)了解別構(gòu)抑制劑對這兩個(gè)子域構(gòu)象的影響和作用機(jī)制，為便于觀察，圖5B放大了這些區(qū)域?？梢钥闯?，對于手指子域，APOC(藍(lán)線)的B因子普遍較小，APOO(黑線)的普遍較大；EBR(紅線)大部分與APOC相當(dāng)甚至偏小，只有小部分(60－70，135－145)較大，與APOO相當(dāng)。這說明無EFV結(jié)合時(shí)，閉合狀態(tài)(APOC)下手指子域的柔性明顯較小，而張開狀態(tài)(APOO)下手指子域總體上表現(xiàn)出較大柔性。對于結(jié)合了EFV的體系EBR，其拇指子域的B因子普遍變小，表明別構(gòu)抑制劑對拇指子域的柔性具有明顯抑制作用，這正是所謂“拇指關(guān)節(jié)炎”效應(yīng)；值得注意的是，這時(shí)手指子域大部分殘基(1－60，70－86)的柔性也變小了，說明EFV對這些殘基也具有抑制作用，手指也患了“關(guān)節(jié)炎”，只是癥狀較輕。對于手掌子域，三個(gè)體系的B因子差別不大，說明手掌子域具有較強(qiáng)的剛性，手指拇指的張合以及EFV的加入對手掌子域骨架支撐結(jié)構(gòu)的影響不大。

圖5　三個(gè)RT體系全部殘基(A)和聚合酶活性區(qū)殘基(B)的B因子Fig.5　B-factors of all the residues of the three RT systems(A)and the polymerase domain(B)

綜上所述，別構(gòu)抑制劑EFV對手掌子域骨架殘基的剛性影響很?。粚δ粗缸佑蚪^大部分殘基的柔性有明顯抑制；對RT手指子域大部分殘基的柔性也有抑制作用。這個(gè)結(jié)果有利地證明了關(guān)節(jié)炎模型。而且，我們的結(jié)果表明EFV不但引起“拇指關(guān)節(jié)炎”，也使手指患上了輕度“關(guān)節(jié)炎”，這是之前文獻(xiàn)沒有報(bào)道的。

3.4構(gòu)象變化

蛋白分子中某些原子間的距離變化常常能反映整個(gè)體系構(gòu)象的變化，這種距離分析方法可以給出直觀清楚的構(gòu)象變化信息。為了理解RT上手形構(gòu)象的動(dòng)力學(xué)性質(zhì)，我們計(jì)算了手指子域頂端Lys66和拇指子域頂端Leu289的Cα原子之間的距離L(圖6A)在動(dòng)力學(xué)過程中的變化。這個(gè)距離由于是拇指和手指遠(yuǎn)端間的距離，因此能夠較好地反映手形的張合程度。

圖6B直觀地給出了動(dòng)力學(xué)過程中三個(gè)體系手形區(qū)域平均結(jié)構(gòu)的疊合示意圖。三個(gè)平均結(jié)構(gòu)均取自RMSD值較小且穩(wěn)定的模擬時(shí)段，分別為，EBR：74－88 ns，紅色；APOO：12－32 ns，黃色；APOC：48－62 ns，藍(lán)色。可以看出，三個(gè)體系的拇指子域結(jié)構(gòu)表現(xiàn)出明顯差異：APOC的拇指子域向手指子域方向大幅度靠攏，發(fā)生握合；相對APOO，EBR的拇指子域發(fā)生了明顯的向外張開偏移。然而比較起來，三體系手指子域的構(gòu)象變化很小。這說明別構(gòu)抑制劑或RT初結(jié)構(gòu)對體系動(dòng)力學(xué)結(jié)構(gòu)的影響主要體現(xiàn)在拇指子域，或者說拇指子域的構(gòu)象對RT的性質(zhì)(結(jié)構(gòu)及由此引起的活性)起著重要作用，這個(gè)直觀圖像與實(shí)驗(yàn)上的“拇指關(guān)節(jié)炎模型”相符。

圖6C清楚地表明，無“楔子”的閉合狀態(tài)下(APOC)，L在整個(gè)動(dòng)力學(xué)過程中基本不變，保持在1.49 nm附近；EBR和APOO兩個(gè)體系的L值大小及波動(dòng)幅度相似，其值都在3.8 nm左右，變化幅度約為±0.5 nm。這個(gè)結(jié)果給出的物理圖像是：在整個(gè)動(dòng)力學(xué)模擬過程中，APOC保持握緊的手形；EBR和APOO則保持張開的手形，而且二者張開的程度相當(dāng)。圖6D給出了三個(gè)體系中L值的幾率(時(shí)間)分布函數(shù)P，它不但反映手形的張合程度，更重要的是它反映了系統(tǒng)構(gòu)象的穩(wěn)定程度。因?yàn)?，一個(gè)確定的L值對應(yīng)一個(gè)確定的構(gòu)象或手形，所以P函數(shù)曲線的形狀對應(yīng)構(gòu)象的變化程度或者結(jié)構(gòu)穩(wěn)定程度。圖6D中，APOC的P(L)分布(藍(lán)線)是一高聳的中心位于1.49 nm的單峰，表明該體系在整個(gè)動(dòng)力學(xué)模擬過程中，處于一種穩(wěn)定的構(gòu)象結(jié)構(gòu)，從自由能曲面上看，代表閉合狀態(tài)對應(yīng)一個(gè)較深的能量極小值點(diǎn)，體系的相點(diǎn)被牢牢的束縛在極值點(diǎn)附近做小幅振動(dòng)。特別是，1.49 nm處的峰不但高聳尖銳而且對稱性極高，這意味無別構(gòu)分子加入時(shí)，握合的手形是一個(gè)非常穩(wěn)定的狀態(tài)，除非有外力作用，在等溫平衡熱力學(xué)條件下不會(huì)自動(dòng)轉(zhuǎn)變成張開狀態(tài)。這就是說，正常情況下，RT總是傾向于閉合的活性狀態(tài)，具有加工合成病毒的功能。這個(gè)結(jié)果與實(shí)驗(yàn)一致，但不同于Ivetac和McCammon15先前的模擬結(jié)果，他們通過改變體系中原子的初始速度使四個(gè)副本中的兩個(gè)發(fā)生了由閉合到張開的狀態(tài)躍遷。值得注意的是，這種構(gòu)象變化都發(fā)生在模擬運(yùn)行的初期，因此可以推斷，這種躍遷與原子初速度的賦值有關(guān)，即這些發(fā)生躍遷的樣本體系相點(diǎn)的初速度剛好處于通向勢壘較低的勢能峽谷中，從而在程序?qū)w系增溫(對原子速度乘以大于1的系數(shù))至300 K時(shí)，推動(dòng)體系越過勢壘進(jìn)入另一能量極值點(diǎn)。另外，計(jì)算誤差和力場誤差也都可能引起不準(zhǔn)確的結(jié)果。EBR體系的P(L)也是一個(gè)對稱性良好，中心位于3.84 nm處的單峰，表明體系也始終穩(wěn)定在相應(yīng)的自由能極值點(diǎn)附近，這正是所謂的“分子楔”作用。比較而言，APOO的P(L)是一個(gè)不對稱分布，在主峰3.74 nm左側(cè)相距約0.30－0.70 nm處出現(xiàn)了一臺梯。從構(gòu)象上說，這個(gè)臺梯結(jié)構(gòu)反映了體系的構(gòu)象由張開向閉合轉(zhuǎn)變的趨勢或征兆。換言之，盡管APOO大部分時(shí)間處于張開的初構(gòu)象，但與EBR不同的是，它存在由張開向閉合轉(zhuǎn)變的明顯傾向，并有相當(dāng)?shù)膸茁侍幱诎胛蘸蠣顟B(tài)，這個(gè)結(jié)果也與實(shí)驗(yàn)一致，即自然的APOO總是趨于具有活性的閉合構(gòu)象。

圖6　100 ns動(dòng)力學(xué)模擬過程中RT手形握合狀態(tài)的變化Fig.6　Conformation changes of the RT hand opening and closing in the 100 ns dynamics simulations

綜上所述，可以看出，相對于RMSD和B因子分析，手指頂端Lys66和拇指頂端Leu289的Cα原子間距離L的分析給出了體系在動(dòng)力學(xué)過程中更加直觀的物理圖像和更加清晰的動(dòng)力學(xué)性質(zhì)。我們的結(jié)果支持基于實(shí)驗(yàn)晶體資料基礎(chǔ)上的“拇指關(guān)節(jié)炎”或“分子楔”模型，并與文獻(xiàn)11,12報(bào)道的模擬結(jié)果一致。第一次揭示了APOO由張開狀態(tài)向閉合構(gòu)象轉(zhuǎn)變的傾向和動(dòng)力學(xué)性質(zhì)，即APOO的張開構(gòu)象也是一個(gè)相對穩(wěn)定的能量極點(diǎn)，只是相應(yīng)于閉合方向的勢能面比較低平，這成為APOO由張向合轉(zhuǎn)變的驅(qū)動(dòng)力。我們沒有得到文獻(xiàn)16報(bào)道的構(gòu)象躍遷現(xiàn)象，原因可能是所用力場不同，也可能是他們使用較大時(shí)間步長積累的誤差所致，更可能與模擬開始時(shí)原子的初速度隨機(jī)賦值及隨后的加溫過程有關(guān)，因?yàn)槲墨I(xiàn)報(bào)道的構(gòu)型轉(zhuǎn)變均發(fā)生在模擬開始階段。

4　結(jié)論

運(yùn)用GPU大規(guī)模并行計(jì)算技術(shù)，采用AMBER12和ff12SB力場，對三種RT大分子體系分別進(jìn)行了100 ns較長時(shí)間的分子動(dòng)力學(xué)模擬，研究了別構(gòu)抑制劑的結(jié)合對體系結(jié)構(gòu)和不同子域殘基柔性的影響。主要結(jié)論是：(1)當(dāng)別構(gòu)抑制劑EFV結(jié)合到RT上的特定別構(gòu)位后，會(huì)導(dǎo)致體系結(jié)構(gòu)的變化，這個(gè)變化會(huì)影響RT的活性，從而實(shí)現(xiàn)別構(gòu)抑制劑對RT功能的控制。(2)通過系統(tǒng)地分析相關(guān)體系的RMSD以及殘基的B因子，發(fā)現(xiàn)對于RT，別構(gòu)抑制劑的作用不但撐開拇指和手指子域，而且增強(qiáng)了兩個(gè)子域殘基的剛性，證實(shí)了“分子楔”模型觀點(diǎn)。值得說明的是我們的結(jié)果不但驗(yàn)證了“拇指關(guān)節(jié)炎”觀點(diǎn)，而且預(yù)測了輕度“手指關(guān)節(jié)炎”癥狀。(3)通過手指頂端Lys66和拇指頂端Leu289的Cα原子間距離L的分析發(fā)現(xiàn)，無別構(gòu)抑制劑加入時(shí)，握合的手形是一個(gè)非常穩(wěn)定的狀態(tài)，即自然條件下，RT總是傾向于閉合的活性狀態(tài)，具有加工合成病毒的功能。EBR體系的張開狀態(tài)也是一個(gè)穩(wěn)定狀態(tài)。比較而言，盡管APOO大部分時(shí)間處于張開的初構(gòu)象，但與EBR不同的是，它存在由張開向閉合轉(zhuǎn)變的明顯傾向，并有相當(dāng)?shù)膸茁侍幱诎胛蘸蠣顟B(tài)，說明自然的APOO傾向于具有活性的閉合構(gòu)象。這些結(jié)果與實(shí)驗(yàn)一致，并證實(shí)了別構(gòu)抑制劑NNRTIs的“分子楔”作用。另外，通過比較不同文獻(xiàn)報(bào)告的模擬結(jié)果，我們認(rèn)為，動(dòng)力學(xué)模擬結(jié)果可能對模擬方法和策略具有很強(qiáng)的依賴性，包括模擬時(shí)間長短、時(shí)間步長的大小、模擬開始結(jié)構(gòu)的準(zhǔn)備、力場的精度以及原子初速度的隨機(jī)賦值(原子的初速度應(yīng)具有熱力學(xué)平衡分布的特點(diǎn)，即Maxwell分布)等，精確的力場和模擬過程中體系的良好動(dòng)力學(xué)平衡是得到可靠模擬結(jié)論的必要條件。最后強(qiáng)調(diào)的是，我們的計(jì)算結(jié)果表明，一旦RT體系穩(wěn)定在開放或閉合狀態(tài)，則狀態(tài)改變是一種小概率事件，常溫下(300 K)僅靠熱力學(xué)效應(yīng)似乎很難使體系在兩態(tài)之間發(fā)生躍遷。對這些問題的闡述將有助于更加詳細(xì)地理解NNRTIs的抑制機(jī)制和聚合酶結(jié)構(gòu)域構(gòu)象變化的動(dòng)力學(xué)性質(zhì)，從而為設(shè)計(jì)更有效的抑制劑提供理論依據(jù)。然而，要得到這些問題的確切結(jié)論，需要進(jìn)行多樣本和更長時(shí)間的分子動(dòng)力學(xué)模擬計(jì)算。希望我們的工作能給其他研究者，尤其是從事模擬計(jì)算的人員帶來啟發(fā)。

References

(1)Mathers,C.D.;Loncar,D.PLoS Med.2006,3(11),e442. (2)De Clercq,E.Chem.Biodivers.2004,1,44.

(3)Ren,J.;Stammers,D.K.Trends Pharmacol.Sci.2005,26, 4.doi:10.1016/j.tips.2004.11.003

(4)Zhu,R.X.;Wang,F.;Liu,Q.;Kang,T.G.Acta Chim.Sin. 2011,69(15),1731.[朱瑞新,王飛,劉琦,康廷國,化學(xué)學(xué)報(bào),2011,69(15),1731.]

(5)Jacobo-Molina,A.;Arnold,E.Biochemistry 1991,30(26), 6351.doi:10.1021/bi00240a001

(6)Lawtrakul,L.;Beyer,A.;Hannongbua,S.;Wolschann,P. Monatsh.Chem.2004,135(8),1033.

(7)Sluis-Cremer,N.;Temiz,N.A.;Bahar,I.Curr.HIV Res.2004, 2(4),323.doi:10.2174/1570162043351093

(8)Bakan,A.;Bahar,I.Proc.Natl.Acad.Sci.U.S.A.2009,106 (34),14349.doi:10.1073/pnas.0904214106

(9)Kohlstaedt,L.A.;Wang,J.;Friedman,J.M.;Rice,P.A.; Steitz,T.A.Science 1992,256(6),1783.doi:10.1126/ science.1377403

(10)Liu,S.X.;Abbondanzieri,E.A.;Rausch,J.W.;Le Grice,S.F. J.;Zhuang,X.W.Science 2008,322(5904),1092.doi: 10.1126/science.1163108

(11)Shen,L.L.;Shen,J.H.;Luo,X.M.;Cheng,F.;Xu,Y.C.; Chen,K.X.;Arnold,E.;Ding,J.P.;Jiang,H.L.Biophys.J. 2003,84(6),3547.doi:10.1016/S0006-3495(03)75088-7

(12)Zhou,Z.;Madrid,M.;Evanseck,J.D.J.Am.Chem.Soc.2005, 127(49),17253.doi:10.1021/ja053973d

(13)Madrid,M.;Jacobo-Molina,A.;Ding,J.;Arnold,E.Proteins 1999,35(3),332.

(14)Madrid,M.;Lukin,J.A.;Madura,J.D.;Ding,J.;Arnold,E. Proteins 2001,45(3),176.

(15)Ivetac,A.;McCammon,A.J.J.Mol.Biol.2009,388(3), 644.doi:10.1016/j.jmb.2009.03.037

(16)Wright,D.W.;Sadiq,S.K.;De Fabritiis,G.;Coveney,P.V. J.Am.Chem.Soc.2012,134(31),12885.doi:10.1021/ ja301565k

(17)Chen,J.Z.;Wang,J.N.;Zhu,W.L.;Li,G.H.J.Comput. Aided Mol.Des.2013,27(11),965.

(18)Luo,F.;Gao,J.;Cheng,Y.H.;Cui,W.;Ji,M.J.Acta Phys.-Chim.Sin.2012,28(9),2191.[羅芳,高劍,成元華,崔巍,計(jì)明娟.物理化學(xué)學(xué)報(bào),2012,28(9),2191.] doi:10.3866/PKU.WHXB201207063

(19)Zhang,H.;Lu,J.R.;Mu,J.B.;Liu,J.B.;Yang,X.Y.;Wang, M.J.;Zhang,R.B.Acta Phys.-Chim.Sin.2015,31(10),566. [張賀,盧俊瑞,穆江蓓,劉金彪,楊旭云,王美君,張瑞波,物理化學(xué)學(xué)報(bào),2015,31(10),566.]doi:10.3866/PKU. WHXB201501061

(20)Case,D.A.;Darden,T.A.;Cheatham,T.E.,III;Simmerling, C.L.;Wang,J.;Duke,R.E.;Luo,R.;Walker,R.C.;Zhang, W.;Merz,K.M.;Roberts,B.P.;Hayik,S.;Roitberg,A.E.; Seabra,G.;Swails,J.M.;Kolossváry,I.;Wong,K.F.;Paesani, F.;Vanicek,J.;Wolf,R.M.;Liu,J.;Wu,X.;Brozell,S.R.; Steinbrecher,T.;Gohlke,H.;Cai,Q.;Ye,X.;Wang,J.;Hsieh, M.J.;Cui,G.;Roe,D.R.;Mathews,D.H.;Seetin,M.G.; Salomon-Ferrer,R.;Sagui,C.;Babin,V.;Luchko,T.;Gusarov, S.;Kovalenko,A.;Kollman,P.A.AMBER 12;University of California:San Francisco,CA,2012.

(21)Meng,X.M.;Wang,J.L.;Zhang,S.L.;Zhang,Q.G.ActaChim.Sin.2013,71(8),1167.[孟現(xiàn)美,王加磊,張少龍,張慶剛,化學(xué)學(xué)報(bào),2013,71(8),1167.]doi:10.6023/A13030327

(22)Lindorff,L.K.;Piana,S.;Palmo,K.;Maragakis,P.;Klepeis,J. L.;Dror,R.O.;Shaw,D.E.Protein Force Field 2010,78(8), 1950.

(23)Chen,J.Z.;Zhang,D.L.;Zhang,Y.X.;Li,G.H.Int.J.Mol. Sci.2012,13(2),2176.

(24)Humphrey,W.;Dalke,A.;Schulten,K.J.Mol.Graphics 1996, 14(1),33.doi:10.1016/0263-7855(96)00018-5

(25)The Theoretical and Computational Biophysics Group.VMD, Revision 1.8.7;NIH Center for Macromolecular Modeling and Bioinformatics,the Beckman Institute,University of Illinois at Urbana-Champaign.

(26)Beyer,A.;Lawtrakul,L.;Hannongbua,S.;Wolschann,P. Monatsh.Chem.2004,135(7),1047.

(27)De Clercq,E.Nat.Rev.Drug Discov.2007,6(12),1001. doi:10.1038/nrd2424

Effect of the Allosteric Inhibitor Efavirenz on HIV-1 Reverse Transcriptase by Molecular Dynamics Simulation

MENG Xian-MeiZHANG Shao-Long*ZHANG Qing-Gang*
(School of Physics and Electronics,Shandong Normal University,Jinan 250014,P.R.China)

To understand the allosteric modulation dynamics of non-nucleoside reverse transcriptase inhibitors (NNRTIs),various models and suggestions have been derived from crystallography and simulation.Here,using a new force field,ff12SB,and GPU parallel computing technology,we performed 100-ns-long molecular dynamics simulations on three reverse transcriptase(RT)systems,one bound to inhibitor Efavirenz(EFV)and the others free.Analyses of the influence of the EFV on the conformation of the RT,flexibility of residues and dynamic behaviors of the systems were conducted.The simulations indicate that EFV binding induces structural distortion of the RT,whereas the configuration of the RT is more stable during dynamics,along with a decreasing extent of motion of the residues.EFV suppresses the flexibility of the thumb subunit and reduces that of most residues in the fingers subdomain as well,suggesting that EFV causes not only the so-called“thumb arthritis”but also a slight“fingers arthritis”.No conformational transition occurred throughout the entire simulations and the samples maintained their starting conformations,i.e.,free RT with a closed conformation stayed in the functional state and EFV-bound RT remained in open conformation.However,EFV-free RT with an initially open conformation exhibited an evident trend toward the closed state.These results agree with the models from experiments,and present a useful insight into the allosteric inhibition mechanism of NNRTIs.In addition,the simulation methodology has been discussed in detail and will be of significance to the computational simulation of large biological molecules.

HIV-1 reverse transcriptase;Nonnucleoside reverse transcriptase inhibitor; Allosteric inhibitor;Molecular dynamics simulation;Conformation

September 14,2015;Revised:November 25,2015;Published on Web:November 30,2015.*Corresponding authors.ZHANG Qing-Gang,Email:zhangqg@sdnu.edu.cn;Tel:+86-531-86180349. ZHANG Shao-Long,Email:slzhang@sdnu.edu.cn;Tel:+86-531-86182521.

O641

The project was supported by the National Natural Science Foundation of China(11274206).

國家自然科學(xué)基金(11274206)資助項(xiàng)目

?Editorial office ofActa Physico-Chimica Sinica