抑制劑8Q9和8QC與bromodomain結(jié)構(gòu)域蛋白4作用機理的分子動力學(xué)研究

吳世亮，孫海波，尹妍妍，李洪云，王 青，王立飛

(山東交通學(xué)院理學(xué)院，濟南 250357)

1 引 言

組蛋白乙酰化主要由所謂的“書寫者”，“閱讀者”和“擦除者”進行調(diào)節(jié)，bromodomain結(jié)構(gòu)域(Bromodomain，BRD)是組蛋白乙?；摹伴喿x者”[1].迄今為止，人類已識別出61種BRDS蛋白結(jié)構(gòu)域，由約140個氨基酸殘基組成.這些結(jié)構(gòu)域存在于46種不同的蛋白質(zhì)中，并分為8個不同的家族[2，3].BRD具有高度保守的折疊模式，主要由四個反向平行的α-螺旋(αA，αB，αC，αZ)和兩個柔性環(huán)(ZA環(huán)和BC環(huán))構(gòu)成[4，5].Bromodomain結(jié)構(gòu)域和末端結(jié)構(gòu)域(bromodomain and extra-terminal domain，BET)家族是八個主要BRD家族的重要組成部分，包含BRD2，BRD3，BRD4和BRDT[6]四個成員.

BRD4是BET家族中重要的功能蛋白.研究表明，BRD4通過招募正向轉(zhuǎn)錄延長因子b(positive transcriptional elongation factor b，P-TEFb)調(diào)節(jié)靶基因的表達，在調(diào)節(jié)細胞基因轉(zhuǎn)錄、細胞周期、炎癥等生物過程中發(fā)揮重要作用[7-10].實際上，肺癌、乳腺癌和肝癌等各種惡性腫瘤的發(fā)生和發(fā)展與BRD4的表達失調(diào)或功能異常密切相關(guān)[11-13].此外，BRD4對急性骨髓性白血病細胞的分化誘導(dǎo)產(chǎn)生重要影響[14].因此，BRD4已成為用于臨床治療各種惡性腫瘤、急性骨髓性白血病、炎性疾病等藥物設(shè)計的重要靶點之一.

自2010年首次公開了有效的BET家族抑制劑JQ1以來[4]，許多針對BET家族特別是靶向BRD4的小分子抑制劑被廣泛報道，如IBET151[15]，I-BET762[16]和RVX-208[17-18]等，這些抑制劑可以誘導(dǎo)細胞周期停滯和凋亡，并抑制腫瘤的侵襲和轉(zhuǎn)移.然而，抑制劑與BRD4蛋白的結(jié)合方式以及蛋白構(gòu)象變化的細節(jié)仍然不足，這限制了針對BRD4的有效抑制劑的設(shè)計與開發(fā).因此，進一步從原子層次了解抑制劑和BRD4的結(jié)合機制變得至關(guān)重要.

隨著計算方法的發(fā)展，分子動力學(xué)模擬[19，20]和結(jié)合自由能計算[21，22]已成為探討抑制劑與蛋白質(zhì)結(jié)合機制的有效工具.在這項工作中，選擇了兩種抑制劑8Q9和8QC來闡明抑制劑與BRD4(1)的結(jié)合機制[23]，我們之所以選擇這兩個小分子抑制劑是因為它們顯示出與BRD4(1)不同的結(jié)合能力，對應(yīng)于8Q9和8QC的IC50值分別為2270和150 nM.圖1展示了抑制劑8Q9、8QC以及抑制劑/BRD4(1)復(fù)合體的結(jié)構(gòu).抑制劑8Q9和8QC與BRD4(1)的結(jié)合差異的闡明對于理解抑制劑與BRD4的結(jié)合機理具有重要意義.因此，本工作采用分子動力學(xué)模擬，結(jié)合自由能計算和抑制劑-殘基相互作用來研究兩種抑制劑與BRD4(1)的結(jié)合方式以及抑制劑結(jié)合導(dǎo)致的BRD4(1)結(jié)構(gòu)的變化.同時，我們也希望這項研究能夠為設(shè)計有效抑制BRD4(1)活性的抑制劑提供有價值的信息和理論指導(dǎo).

圖1 分子的結(jié)構(gòu)：(a)抑制劑8Q9與BRD4(1)復(fù)合物的結(jié)構(gòu)，(b)抑制劑8Q9的結(jié)構(gòu)，(c)抑制劑8QC的結(jié)構(gòu).Fig.1 Structures of molecules：(a)the8Q9 and BRD4(1)complex，(b)the inhibitor 8Q9，(c)the inhibitor 8QC.

2 理論和方法

2.1 模擬系統(tǒng)的準備

動力學(xué)模擬的初始模型取自蛋白質(zhì)庫中的晶體結(jié)構(gòu)5Y93和5Y94[23]，5Y93是抑制劑8Q9和BRD4(1)的復(fù)合物(8Q9/BRD4(1))，5Y94是抑制劑8QC和BRD4(1)的復(fù)合物(8QC/BRD4(1))，apo/BRD4(1)表示從晶體結(jié)構(gòu)5Y94中去掉抑制劑8QC獲得的BRD4(1)的無抑制劑狀態(tài)，該結(jié)構(gòu)用于比較抑制劑結(jié)合對BRD4(1)構(gòu)象變化的影響.所有結(jié)晶水分子都保留在初始模型中，借助Amber18中的LEAP模塊，將晶體結(jié)構(gòu)中缺失的氫原子連接到相應(yīng)的重原子.對兩種小分子抑制劑，使用了Amber18的antechamber程序中的AM1方法來實現(xiàn)結(jié)構(gòu)優(yōu)化和計算，并將原子的BCC電荷分配給抑制劑的各原子.蛋白質(zhì)和水分子的力場參數(shù)由Amber18中的ff14SB力場產(chǎn)生，而小分子抑制劑的力場參數(shù)由GAFF力場產(chǎn)生.復(fù)合體被溶解在顯性的TIP3P水盒子里，水盒子邊緣與復(fù)合體最近原子的距離設(shè)定為12.0?，然后將三個氯離子置于水和復(fù)合體組成的系統(tǒng)中，以確保系統(tǒng)呈電中性.

2.2 分子動力學(xué)模擬

對于每個復(fù)合體系統(tǒng)，通過Amber18中的PMEMD模塊執(zhí)行了能量最小化和分子動力學(xué)模擬.為了消除一些原子間的不合理接觸，每個系統(tǒng)都經(jīng)過了兩階段的能量優(yōu)化.首先，對溶質(zhì)進行約束并優(yōu)化溶劑.其次，整個系統(tǒng)實現(xiàn)了無約束的優(yōu)化.在優(yōu)化過程中，首先執(zhí)行2500步的最陡下降優(yōu)化，接著進行2500步的共軛梯度優(yōu)化.之后，將每個系統(tǒng)在2 ns內(nèi)保持恒定體積逐漸從0 K加熱到300 K，接著在常溫300 K，常壓1標準大氣壓條件下進行1 ns分子動力學(xué)平衡.最后，進行200 ns的無約束分子動力學(xué)模擬.在模擬過程中，每2 ps記錄一次原子坐標.使用SHAKE算法來限制所有涉及氫原子的化學(xué)鍵的伸縮，模擬積分步長為2 fs.借助Langevin恒溫器控制每個系統(tǒng)的溫度.PME方法用于計算原子之間的長程靜電相互作用.涉及靜電相互作用和范德華相互作用的非成鍵相互作用計算的截斷值設(shè)置為10.0?.

2.3 結(jié)合自由能的計算

MM-GBSA方法是計算結(jié)合自由能的一種常用方法.在我們的研究中，借助Amber18中的MM-GBSA方法有效地估計了抑制劑與BRD4(1)的結(jié)合自由能.對于每個系統(tǒng)，從動力學(xué)模擬軌跡的最后80 ns中每隔400 ps取一個構(gòu)象，共取出200個構(gòu)象用于MM-GBSA方法的計算.計算過程中去掉構(gòu)象中的水分子和離子，采用下面的公式計算抑制劑和蛋白質(zhì)的結(jié)合自由能：

式中△Eele和△Evdw分別表示氣相中抑制劑和蛋白質(zhì)的靜電相互作用與范德華相互作用，這兩項用分子力學(xué)計算，力場參數(shù)由Amber18中的ff14SB力場產(chǎn)生；第三項△Gpol表示極性溶解自由能對抑制劑和蛋白質(zhì)結(jié)合的貢獻，通常由Amber套件中MM-GBSA模塊求得；第四項△Gnonpol表示非極性溶解自由能對抑制劑和蛋白質(zhì)結(jié)合的貢獻，使用下面的方程求解：

式中SASA是溶劑可及表面積，經(jīng)驗常數(shù)γ和β的值分別設(shè)置為0.0072 kcal·mol-1·?-2和0.0 kcal·mol-1；最后一項-T△S表示熵變對抑制劑和蛋白質(zhì)結(jié)合的貢獻，是平動，轉(zhuǎn)動和振動熵之和，該項使用振動模和傳統(tǒng)的熱力學(xué)計算.

2.4 動態(tài)相關(guān)性分析

為了研究抑制劑結(jié)合引起的BRD4(1)的運動模式的改變，采用Amber 18中的CPPTRAJ程序，基于分子動力學(xué)模擬的平衡部分，對三個系統(tǒng)進行動態(tài)相關(guān)性分析.由第i個Cα原子和第j個Cα原子構(gòu)成的矩陣元素Cij定義如下：

其中△ri表示第i個Cα原子相對于其平均值的位置偏差.Cij的取值范圍是-1至+1，其中-1表示完全負相關(guān)運動，而+1表示完全正相關(guān)運動.正相關(guān)說明殘基沿相同方向移動，而負相關(guān)說明殘基沿相反方向移動.使用顏色編碼模式來體現(xiàn)殘基之間相關(guān)運動的程度.文中出現(xiàn)的DCCM圖是使用分子動力學(xué)模擬軌跡的最后80 ns中保存的Cα原子的坐標計算得到的.

3 結(jié)果和討論

3.1 分子動力學(xué)平衡

為了評估動力學(xué)平衡的穩(wěn)定性，分析了三個系統(tǒng)的BRD4(1)主鏈原子相對于其初始優(yōu)化結(jié)構(gòu)的均方根偏差(RMSD)隨時間的演化情況(圖2).如圖2所示，三個系統(tǒng)的RMSD波動不同，apo/BRD4(1)，8Q9/BRD4(1)，8QC/BRD4(1)達到穩(wěn)定平衡的時間分別是模擬開始后10 ns，50 ns和120 ns，其平衡后的平均RMSD值分別為1.76，0.98和1.18?，三個系統(tǒng)的RMSD的漲落范圍均小于0.65?，該結(jié)果說明分子動力學(xué)平衡的穩(wěn)定性是可靠的，可用于隨后分析的統(tǒng)計取樣.此外，8Q9/BRD4(1)，8QC/BRD4(1)兩個復(fù)合體的RMSD值和漲落范圍明顯低于apo/BRD4(1)系統(tǒng)，表明抑制劑結(jié)合對BRD4(1)的結(jié)構(gòu)性波動產(chǎn)生某些限制，這有助于抑制劑/BRD4(1)復(fù)合體的動力學(xué)穩(wěn)定性.

圖2 分子動力學(xué)模擬中BRD4(1)主鏈原子的均方根偏差Fig.2 Root-mean-square deviations(RMSDs)of the backbone atoms in BRD4(1)relative to the corresponding crystal structures as function of simulations time.

3.2 結(jié)構(gòu)變化

殘基的Cα原子均方根漲落(Root-meansquare fluctuations，RMSFs)可用于評估抑制劑結(jié)合對BRD4(1)結(jié)構(gòu)柔性的影響.采用Amber18中的CPPTRAJ工具計算了apo/BRD4(1)，8Q9/BRD4(1)，8QC/BRD4(1)三個系統(tǒng)中各殘基Cα原子漲落(圖3).圖3表明抑制劑的存在對BRD4(1)的結(jié)構(gòu)柔性產(chǎn)生了重要影響.抑制劑結(jié)合的BRD4(1)中殘基89-97和136-145兩個區(qū)域的RMSF值比apo/BRD4(1)的要小，這表明這兩個區(qū)域的柔性由于抑制劑的結(jié)合而減弱.這些結(jié)果表明，上述區(qū)域中的一些關(guān)鍵殘基可能與兩種抑制劑發(fā)生強烈的相互作用，這可以通過后面基于殘基的自由能分解和氫鍵分析的計算進一步證明.

圖3 動力學(xué)模擬過程中三個系統(tǒng)的BRD4(1)殘基Cα原子的平均漲落Fig.3 Root-mean-square fluctuations(RMSFs)of Cα atoms in BRD4(1)of three simulated systems

3.3 動態(tài)相關(guān)性分析

為了研究抑制劑結(jié)合引起的BRD4(1)的運動模式的改變，采用Amber 18中的CPPTRAJ程序，基于分子動力學(xué)模擬的平衡部分，計算了殘基之間的DCCM(圖4)，圖中借助顏色編碼模式來體現(xiàn)特定殘基之間相關(guān)運動的程度，對角區(qū)域?qū)?yīng)于殘基相對于它們自身的運動，而非對角區(qū)域代表不同殘基之間的相對運動，紅色和黃色表示殘基之間的強烈正相關(guān)運動，而藍色和深藍色表示殘基之間的強烈負相關(guān)運動.

apo/BRD4(1)的DCCM如圖4(a)所示，從圖中可以看出在區(qū)域R1和R3(黃色)中發(fā)生了強烈的正相關(guān)運動，區(qū)域R1中的運動指示殘基126-140相對于殘基99-108的相關(guān)運動，而在R3區(qū)域中的運動描述殘基116-121和54-59之間的相關(guān)運動；在區(qū)域R2(綠色)中注意到了一些輕微的正相關(guān)運動，描述了殘基81-89相對于殘基102-110的弱相關(guān)運動.與apo/BRD4(1)相比，8Q9/BRD4(1)中抑制劑的存在增強了R1和R2區(qū)域中發(fā)生的正相關(guān)運動，而大大減弱了R3區(qū)域的正相關(guān)運動，此外，8Q9的結(jié)合明顯增強了區(qū)域R4中的正相關(guān)運動(圖4(b)).圖4(c)表明8QC/BRD4(1)中抑制劑的存在減弱了R1和R2區(qū)域中發(fā)生的正相關(guān)運動，而增強了R3區(qū)域的正相關(guān)運動，此外，8QC的結(jié)合明顯減弱了區(qū)域R5中的負相關(guān)運動.以上分析表明，抑制劑的結(jié)合會導(dǎo)致BRD4(1)的運動模式發(fā)生重大變化，進而可能會對抑制劑與BRD4(1)中各個分離殘基的相互作用產(chǎn)生重大影響.

圖4 動態(tài)相關(guān)性分析：(a)apo/BRD4(1)，(b)8Q9/BRD4(1)，(c)8QC/BRD4(1).Fig.4 Dynamic cross-correlation maps(DCCMs)：(a)apo/BRD4(1)，(b)8Q9/BRD4(1)，(c)8QC/BRD4(1).

3.4 結(jié)合自由能和作用機理

為揭示影響抑制劑結(jié)合的主要因素，比較了兩種抑制劑和BRD4(1)的結(jié)合親和力的差異.在我們的研究中，通過使用MM-GBSA方法有效地估計了兩種抑制劑對BRD4(1)的結(jié)合自由能.對于每個結(jié)合抑制劑的系統(tǒng)，從記錄在最后80 ns的分子動力學(xué)軌跡中以400 ps的時間間隔提取了200個構(gòu)象用于計算.表1中的結(jié)果顯示計算的結(jié)合自由能的趨勢與實驗值的趨勢一致，這表明我們目前的自由能分析是可靠的.如表1所示，結(jié)合自由能分為五個部分，即范德華相互作用(△Evdw)，靜電相互作用(△Eele)，極性溶解自由能(△Gpol)，非極性溶解自由能(△Gnopol)和熵變貢獻(-T△S).

表1 MM-GBSA方法計算的抑制劑與BRD4(1)的結(jié)合自由能Table 1 Binding free energies of inhibitors to BRD4(1)calculated by MM-GBSA methoda

表1表明范德華相互作用，靜電相互作用和非極性溶解自由能有利于抑制劑結(jié)合，而極性溶解自由能和熵變貢獻不利于抑制劑的結(jié)合.抑制劑8Q9和8QC與BRD4(1)的范德華相互作用分別為-40.68±0.26，-37.54±0.24 kcal·

mol-1，這為抑制劑與BRD4(1)的結(jié)合提供了有利的因素.抑制劑8Q9和8QC與BRD4(1)的靜電相互作用分別為-17.34±0.39，-22.32±0.57 kcal·mol-1，8Q9/BRD4(1)和8QC/BRD4(1)的極性溶解自由能分別為35.43±0.42、33.61±0.48 kcal·mol-1，可以看到不利的極性溶解自由能完全屏蔽了有利的靜電相互作用，進一步產(chǎn)生了不利的相互作用(△Gele+pol)，從而嚴重削弱了兩種抑制劑與BRD4(1)的結(jié)合.此外，盡管非極性溶解自由能也為抑制劑結(jié)合提供了有益的作用力，但非極性溶解自由能比范德華相互作用弱得多.因此，范德華相互作用在抑制劑與BRD4(1)的結(jié)合中起著重要作用.

3.5 抑制劑-殘基相互作用

基于殘基的自由能分解的分析可以闡明單個殘基在抑制劑結(jié)合中的作用.為了進一步了解各個分離殘基與抑制劑的結(jié)合，將兩個抑制劑與BRD4(1)復(fù)合體的總結(jié)合自由能分解為每個殘基的貢獻(圖5)，并且采用Amber18中的CPPTRAJ模塊分析了抑制劑與BRD4(1)的氫鍵相互作用(表2).圖6和圖7采用動力學(xué)模擬軌跡生成的最低能級結(jié)構(gòu)，描繪了關(guān)鍵殘基相對于抑制劑的幾何位置.

表2 動力學(xué)模擬過程中抑制劑與關(guān)鍵殘基形成的氫鍵Table 2 The hydrogen bonds formed between key residues and inhibitors.

圖5 抑制劑與BRD4(1)中各分離殘基的相互作用譜：(a)8Q9/BRD4(1)，(b)8QC/BRD4(1)Fig.5 Interactions between two inhibitors and separate residues of BRD4(1)：(a)the 8Q9/BRD4(1)，(b)the 8QC/BRD4(1).

圖6 抑制劑相對于BRD4(1)的關(guān)鍵殘基的幾何位置：(a)8Q9/BRD4(1)，(b)8QC/BRD4(1).Fig.6 Geometric positions of inhibitors relative to key residues of BRD4(1)involving strong interaction：(a)8Q9/BRD4(1)，(b)8QC/BRD4(1).

圖7 抑制劑與BRD4(1)中關(guān)鍵殘基的氫鍵相互作用：(a)8Q9/BRD4(1)，(b)8QC/BRD4(1).Fig.7 Hydrogen bonds formed between two inhibitors and key residues of BRD4(1)：(a)8Q9/BRD4(1)，(b)8QC/BRD4(1).

圖5(a)表明六個殘基與8Q9的相互作用強于1 kcal·mol-1，這些殘基是Trp81，Pro82，Val87，Leu92，Leu94和Ile146.Ile146是這些殘基中和抑制劑相互作用最強的殘基，相互作用能是-3.30 kcal·mol-1，主要來自Ile146烷基與抑制劑的疏水環(huán)的CH-π相互作用(圖6(a)).殘基Trp81和Pro82與抑制劑的相互作用能分別是-1.48和-1.34 kcal·mol-1，該相互作用源自殘基的疏水環(huán)臨近抑制劑的疏水環(huán)從而產(chǎn)生的π-π相互作用.殘基Val87，Leu92和Leu94因其各自的烷基靠近抑制劑的疏水環(huán)從而產(chǎn)生CH-π相互作用，分別貢獻了-1.70，-2.14和-1.57 kcal·mol-1的作用能.

圖5(b)表明六個殘基與8QC的相互作用強于1 kcal·mol-1，這些殘基是Trp81，Pro82，Lys91，Leu92，Cys136和Ile146，其中殘基Ile146和抑制劑相互作用最強，相互作用能是-3.45 kcal·mol-1，該相互作用在結(jié)構(gòu)上吻合8QC的疏水環(huán)與Ile146烷基的CH-π相互作用(圖6(b)).殘基Trp81和Pro82與抑制劑的相互作用能分別是-1.14和-1.59 kcal·mol-1，該相互作用是兩個殘基的疏水環(huán)與抑制劑的疏水環(huán)產(chǎn)生的π-π相互作用.殘基Leu92的烷基靠近抑制劑的疏水環(huán)從而產(chǎn)生CH-π相互作用，且貢獻了-2.23 kcal·mol-1的相互作用能.上述分析表明疏水相互作用(例如CH-π和π-π相互作用)在抑制劑與BRD4(1)的結(jié)合中起關(guān)鍵作用.

氫鍵對于抑制劑與蛋白質(zhì)的結(jié)合是必不可少的.在我們的研究中，基于最后80 ns的分子動力學(xué)模擬軌跡，采用Amber18中的CPPTRAJ模塊，計算了抑制劑與關(guān)鍵殘基間的氫鍵距離，角度和占有率(表2和圖7).如表2所示，抑制劑8Q9與殘基Asn140形成兩個氫鍵相互作用，其對抑制劑結(jié)合的能量貢獻為-0.81 kcal·mol-1(圖5(a))，兩個氫鍵的占有率分別為95.75和56.03%，根據(jù)圖7，該氫鍵作用源自殘基Asn140的NH基團與8Q9的疏水性環(huán)中的氧原子(OAG)及氮原子(NAH)的相互作用；抑制劑8QC與殘基Gln85，Lys91和Cys136形成三個氫鍵相互作用，其對抑制劑結(jié)合的能量貢獻分別為-0.76，-1.62和-1.41 kcal·mol-1(圖5(b))，三個氫鍵的占有率分別為41.71，59.73和56.46%，根據(jù)圖7，三個氫鍵作用分別源自殘基Gln85的NH基團與8QC的氧原子(OAM)，殘基Lys91的氧原子與8QC的NH基團，以及殘基Cys136的SH基團與8QC的氮原子(NAH)的相互作用.

鑒于上述分析，疏水性相互作用(例如π-π相互作用和CH-π相互作用)在抑制劑與BRD4(1)的結(jié)合中起關(guān)鍵作用.此外，氫鍵相互作用也是抑制劑與BRD4(1)結(jié)合的重要因素.常見殘基Trp81，Pro82，Gln85，Val87，Lys91，Leu92，Leu94，Cys136，Asn140和Ile146利于抑制劑與蛋白質(zhì)的結(jié)合，這些殘基在設(shè)計針對BRD4(1)的有效抑制劑時需要特別關(guān)注.因此，在開發(fā)針對BRD4(1)的高效抑制劑時，應(yīng)特別注意疏水相互作用和氫鍵相互作用.

4 結(jié) 論

本研究對apo/BRD4(1)，8Q9/BRD4(1)，8QC/BRD4(1)三個系統(tǒng)進行了200 ns的分子動力學(xué)模擬，以揭示抑制劑與BRD4(1)的結(jié)合方式.RMSFs結(jié)果表明抑制劑的結(jié)合對BRD4(1)的結(jié)構(gòu)柔性產(chǎn)生了重大影響.動態(tài)相關(guān)性分析表明抑制劑的結(jié)合會導(dǎo)致BRD4(1)的運動模式發(fā)生重大變化.借助MM-GBSA方法有效地估算了兩種抑制劑對BRD4(1)的結(jié)合自由能，計算結(jié)果的趨勢與實驗值吻合良好.同時，結(jié)合自由能結(jié)果表明范德華相互作用在抑制劑與BRD4(1)的結(jié)合中起著重要作用.自由能分解的結(jié)果表明殘基Trp81，Pro82，Val87，Leu92，Leu94和Ile146主要與抑制劑產(chǎn)生疏水相互作用(例如π-π相互作用和CH-π相互作用)，而殘基Gln85，Lys91，Cys136和Asn140與抑制劑產(chǎn)生了氫鍵相互作用，因此，這些殘基在今后設(shè)計靶向BRD4(1)抑制劑時應(yīng)特別注意.綜上所述，以上結(jié)果表明控制抑制劑與BRD4(1)結(jié)合的主要是疏水性相互作用和氫鍵相互作用.我們希望這項研究可以為以后開發(fā)抑制BRD4活性的有效抑制劑提供有用的信息.