基于Discovery Studio分子模擬藥物設計軟件的拉伸分子動力學

張心苑，毛雪石

（北京協和醫學院中國醫學科學院藥物研究所，北京 100050）

0 引言

Discovery Studio （簡稱DS）是基于Windows/Linux系統，面向藥學領域的新一代分子建模和模擬環境，它應用于藥物結構功能研究以及藥物發現[1]。Discovery Studio研究領域包括：疾病的發病機理研究、新藥發現和設計、計算藥物化學、腫瘤藥物研究等[2]。

Discovery Studio功能模塊包括基于結構的藥物發現和設計模塊、基于片段的藥物設計模塊、基于藥效團的藥物發現和設計模塊、基于小分子的藥物發現和設計模塊蛋白質建模及模擬模塊、分子力學、分子動力學、分子力學/量子力學模擬模塊等。Discovery Studio 2017新增拉伸分子動力學功能，在分子動力學的基礎上引入拉伸力場，能夠更好地分析配體（藥物分子）、受體的結合方式，本文主要介紹分子動力學模塊的拉伸分子動力學功能。

拉伸分子動力學研究的主要問題是利用拉伸分子動力學計算配體（藥物分子）結合受體的自由能、研究配體（藥物分子）的解離過程、研究蛋白質的去折疊過程或構象變化、拉伸離子或配體（藥物分子）通過通道蛋白、構建平均力勢（PMF）。蛋白質類別（結構）不同，展開不同結構的蛋白質所需的拉伸力、溫度不一樣，拉伸分子動力學過程（由高溫和外力誘導的去折疊過程）有差異[3]。不同蛋白質的拉伸分子動力學過程的差異行為與蛋白質折疊的拓撲結構和蛋白質的二級結構元素的穩定性有關。

1 建立并運行拉伸分子動力學模擬

本文以1TEN、1TIT的晶體結構為例，在Discovery Studio中載入1TEN、1TIT的pdb文件。其中，1TEN是一種纖維粘連蛋白，1TIT是一種免疫球蛋白，兩種蛋白的二級結構都是β折疊、β轉角[4-6]。

1.1 選擇一對原子作為拉伸方向

用rNc表示蛋白質一級結構（氨基酸序列）中N端（NH2-）中氮原子到C端（—COOH）中碳原子連線所形成的向量。拉伸力的方向平行于rNc并且力的方向是使N端（NH2-）中氮原子到C端（—COOH）中碳原子遠離的。

首先，確定N端及C端氨基酸，右鍵單擊molecule window出現display style，選擇protein，再選擇color by,選擇N-to-C Terminal，出現的3D結構視圖顯示只有藍、紅兩部分。在藍色部分內的氨基酸就是N端的氨基酸，在紅色部分內的氨基酸就是C端的氨基酸。N-C端向量即為拉伸方向，蛋白質的結構特點不同確定的拉伸方向不同。

1.2 設置拉伸動力學參數，建立拉伸動力學模擬

1.2.1 確定運算時間

如果恒力足夠大，蛋白質拉伸運算完成的時間與合適的力拉伸運算完成的時間不同[7]。運算時間的確定與施加力的大小有關，需要調節拉伸力的同時觀察拉伸分子動力學運算時間是否合適。從標準分子動力學運算結果可看出1TIT蛋白質中運算時間的平衡點在2400ps后，因此將運算時間設為3000ps最為合適。

1.2.2 確定運算溫度

最佳的拉伸分子動力學溫度為450K，在300K-450K的溫度范圍內都合適，溫度為500K時，拉伸時間小于500ps，拉伸速率較快，不適合觀察拉伸分子動力學過程中蛋白質構象的變化[8]。

1.2.3 確定運算拉伸力

拉伸力等于150pN時可以達到70-250?的拉伸距離，無法完全拉伸蛋白質[9]。拉伸力小于等于250pN時無法完全拉伸蛋白質，拉伸力大于400pN時拉伸速度過快，不適于觀察蛋白質結構的變化，最佳拉伸力為400pN。

1.2.4 運行拉伸分子動力學模擬

運算拉伸動力學后出現600個conformation（構象），蛋白質完全拉伸后的分子構象如圖1所示。

圖1 拉伸方向圖（左） 拉伸分子動力學運算結果圖（右）

2 拉伸動力學結果分析

2.1 觀察不同時間點、不同拉伸距離的蛋白質結構狀態

拉伸動力學方法可以提供整個過程中原子水平的細節，從軌跡中獲得一些快照，可表示拉伸動力學過程中的構象變化?，F分別在rNC=43.26?，rNC=56.99?，rNC=64.02?，rNC=71.72?，rNC=79.96?，rNC=87.15?，rNC=94.83?，rNC=102. 57? 及初始狀態和完全拉伸狀態按照時間節點將8幅圖顯示如圖2。（1 ?=10-10m，?主要用于原子物理、晶體學）。

圖2 不同時間點的軌跡圖

2.2 繪制拉伸距離隨時間變化圖，觀察拉伸動力學的平衡點

觀察拉伸分子動力學模擬結果，從模擬的第一個構象一直到最后一個構象的變化動畫，可以發現距離量是一直變化的，拉伸距離的變化范圍為43.26 ?-102.57 ?。

分析拉伸距離隨時間變化的原因是找到拉伸平衡狀態，最佳狀態時的拉伸距離隨時間變化圖應該是逐漸上升趨勢增加到一定值時趨于平衡的狀態[10]。若增加拉伸模擬時間或增加拉伸力可以更快地達到平衡狀態，做出含平衡狀態的拉伸距離隨時間變化圖。

2.3 繪制RMSD，編譯動力學軌跡監控數據

RMSD圖表示蛋白質結構差異，趨勢逐漸增大至穩定，如圖3所示。Tools explorer選項卡中選擇Analyze Trajectory，再選擇 Analyze，Create Plot設置RMSD參數：Properties選擇Distance1，Analysis type選擇RMSD、RMSF、Properties，其他參數按默認設置即可。

2.4 根據多條運算軌跡計算自由能

根據拉伸分子動力學運算結果適當調整拉伸力、拉伸時間、拉伸溫度的值，會得到若干條運算軌跡，根據30-40條軌跡可以準備地計算自由能，現根據一條拉伸分子動力學計算自由能。Tools explorer選項卡中展開Run Simulations，在Steered Molecular Dynamics菜單下點擊Calculate Free Energy from Steering Force，得到計算自由能隨拉伸距離變化的圖4。

圖3 RMSD圖

圖4 自由能隨拉伸距離的變化圖

3 結論

在Discovery Studio 2017軟件中新增拉伸分子動力學的計算功能，通過軟件模擬拉伸分子動力學過程體現蛋白質的去折疊過程中的構象變化、配體（藥物分子）的結合及解離過程，進一步研究并優化藥物分子與受體蛋白的結合方式。Discovery Studio 2017軟件的拉伸分子動力學模塊有強大的分析動力學運算結果的工具，能夠有效分析藥物分子在力場中與靶標蛋白的結合方式。拉伸分子動力學的計算功能還有利于基于靶標的藥物設計及藥物的高通量篩選，在模擬大分子藥物、化學、生物系統中具有重要作用。

[1] 李亞梅,廖端芳.應用Discovery Studio軟件預測煙酸姜黃素酯抗動脈粥樣硬化的作用靶點[J].中國醫藥導報,2013,10(17):16-18.

[2] Paci E,Karplus M.Unfolding Proteins by External Forces and Temperature:The Importance of Topology and Energetics[J].Proceedings of the National Academy of Sciences of the United States of America,2000,97(12):6521-6526.

[3] 王浩,徐利楠,孫玉娜,等.Ydj1p二聚體中β14～β15與domain-Ⅲ分離的拉伸分子動力學模擬研究[J].生物信息學,2013,11(3):167-171.

[4] Wang Y,Cai W S,Chen L,et al.Molecular dynamics simulation reveals how phosphorylation of tyrosine 26 of phosphoglycerate mutase 1 upregulates glycolysis and promotes tumor growth[J].Oncotarget,2011,8(7).

[5] 季海濤,張萬年,周有駿,等.藥物化學計算機輔助教學方法的建立與應用[J].中國醫學教育技術,2002,16(3):169-170.

[6] Paci E,Caflisch A,Plückthun A,et al.Forces and energetics of hapten-antibody dissociation: a biased molecular dynamics simulation study[J].Journal of Molecular Biology,2001,314(3):589-605.

[7] Boyce,S. E.;Mobley D. L.;Rocklin G.J.;et al. Predicting Ligand Binding Affinity with Alchemical Free EnergyMethodsin a Polar Model Binding Site[J].J.Mol.Biol,2009,394,747-763.

[8] 周雷,徐利楠,周小紅,等.Hsc70-Auxilin復合物的拉伸分子動力學模擬[J].生物信息學, 2014(2):110-116.

[9] Bavi N,Bavi O,Vossoughi M,et al.Nanomechanical properties of MscL alpha helices: A steered molecular dynamics study[J].Channels,2016,11(3).

[10] Park M S.Molecular Dynamics Simulations of the Human Glucose Transporter GLUT1[J].Plos One,2015,10(4):e0125361.