Nadine Homeyer and Holger Gohlke, 原始文檔

2018-03-26 21:44:05 翻譯: 周盛福

請注意: 此處介紹的FEW工具可以在AmberTools 14或更高版本中獲取, 並可與AmberTools和AMBER 14或更高版本一起使用. AmberTools和AMBER 14可以按照AMBER手冊中的安裝指導進行安裝. 在安裝過程中會自動修復BUG並獲取更新補丁.

FEW工具可用於自動配置與同一受體結合的一系列配體的結合自由能計算, 本教程旨在演示該功能. FEW工具支持的設置自由能計算的方法有: 基於分子力學的泊松-玻爾茲曼表面積(MM-PBSA)模型、基於分子力學的廣義波恩表面積(MM-GBSA)模型, 線性相互作用能(LIE)模型和熱力學積分(TI)模型方法. 在使用FEW工具之前, 您應該熟悉在AMBER中運行基本的分子動力學(MD)模擬以及您希望採用的自由能計算方式的理論和設置. 演示分子動力學模擬的設置與運行和通過MM-PB(GB)SA及TI方法計算自由能的教程可以在http://ambermd.org/tutorials中找到. 關於自由能計算方法的理論解釋可以參考相關綜述和書籍: MM-PB(GB)SA [1, 2], LIE [3, 4]和TI [5, 6]. 此外, 強烈建議閱讀文獻N. Homeyer, H. Gohlke, FEW - A Workflow Tool for Free Energy Calculations of Ligand Binding. J. Comput. Chem. 2013, 34, 965-973. 本文中演示的自由能計算示例就是選自這篇文獻介紹的研究案例. 我們確定了3種在凝血級聯反應中起重要作用, 抑制蛋白酶因子Xa的配體的相對結合自由能(表1).

需要的輸入結構:

受體的PDB結構, 即因子Xa蛋白

已將坐標定義到結合口袋的配體結構(mol2格式)

如果配體的結合模式是未知的, 則可以通過類似配體的結合模式構建或通過對接程序獲取. 在這裡顯示的樣本分析中, 使用了前一策略. 請注意, 結合自由能預測的準確度隨著蛋白質結構質量的下降和配體結合模式的不確定性而降低. 因此, 應該使用高質量的蛋白質(即受體)結構, 並且配體的結合位置應儘可能準確地構建.

表1 因子Xa的配體結構和實驗檢測的Ki值(nM)

準備工作

為了便於在本教程中使用FEW工具, 請設置一個變數指定FEW在您系統上的安裝路徑, 例如:

setenvFEW = /home/src/FEW (對於csh或tcsh)或

exportFEW = /home/src/FEW (用於bash、zsh、ksh等)

在定義環境變數後執行, 就應該顯示FEW安裝路徑.

本教程中使用的輸入文件可以在這裡下載. 在解壓文件之後, 按照README文件中的安裝說明, 本教程所需的文件應位於目錄下.

現在, 請創建一個新的文件夾, 並進入到此文件夾中. 在本教程中, 我們假設新文件夾名為, 位於. 以下教程中, 您需要根據您的系統上創建的tutorial文件夾的位置來調整路徑.

現在, 將文件夾(包含所有需求的輸入文件), (包含Mol2格式的配體結構)和(包含FEW的命令文件)複製到您創建的tutorial文件夾中:

cp-r$FEW/examples/tutorial/input_info .

cp-r$FEW/examples/tutorial/structs .

cp-r$FEW/examples/tutorial/cfiles .

現在, 在本地文件夾中我們已經有了所有需要的數據, 因此我們可以開始使用FEW工具來設置自由結合能計算了.

第一節: 分子動力學(MD)模擬的設置

MM-PB(GB)SA和LIE計算是基於MD模擬的軌跡快照進行的. 因此, 在開始這些類型的自由能計算之前, 需要執行該步驟. TI計算可以直接從晶體結構開始, 但通常建議首先運行一個短的MD平衡來消除結構中不好的接觸. 因此, 只有後面的步驟將在這裡演示. 對於之前的步驟, 請參閱AMBER手冊的FEW部分.

在這裡, MD模擬被分成兩步設置, 首先確定配體的參數(步驟1A), 然後產生MD模擬的輸入文件(步驟1B). 圖1展示了這兩個步驟所需的輸入信息和FEW生成的輸出信息.

注意: 為了演示的目的, 在本教程中將分別設置複合物、受體和配體的MD模擬(3-軌跡方法). 如果您希望根據「1-軌跡方法」只進行MM-PB(GB)SA計算(即無LIE或TI計算), 則建議此時只設置複合物的模擬.

圖1: 使用FEW進行MD模擬設置的示意圖. 本教程中使用的輸入信息以藍色顯示, 對於其他體系/參數設置時, 可能需要的輸入數據以灰色顯示, FEW的命令文件以紅色字體標出. 本例中未使用到的輸入文件(上圖中的灰色部分)可以在文件夾``$FEW/examples/input_info中找到, 並且AMBER手冊的FEW部分提供了這些數據的解釋.

A: 參數文件的準備

MD模擬的設置是從定義所需的參數以及重新格式化輸入結構文件開始的. 為此, 我們使用文件夾中的命令文件.

文件可以直接用作FEW的輸入命令文件. 只有紅色的基本輸入/輸出路徑需要根據您系統上的tutorial文件夾的位置進行調整. 為了清晰起見, 命令文件中的注釋以綠色顯示, 只有關鍵設置項和參數以黑色字體給出.

在命令文件中, 第一部分參數定義了基本輸入和輸出數據. 由於在本教程中, 我們使用單個的MOL2結構, 因此不需要設置multi_structure_ligand_file項. bound_rec_structure也沒有指定, 因為這裡的設置中, 在有配體結合的複合物和自由的受體兩種情況下, 我們只使用了一個蛋白質結構. 如果受體結構在配體結合後構象發生了改變, 則可以提供兩種不同的受體結構, 一個是結合狀態(bound_rec_structure)、一個是自由狀態(rec_structure). 我們在此使用的受體結構來源於蛋白質資料庫中的晶體結構2RA0. 由於保留了提交的受體結構(2RA0_IN.pdb)中鑒定的結合水分子(定義為HOH), 所以water_in_rec項設為1.

命令文件中的第二部分包含生成LEaP基本輸入文件所需的參數. 這裡需要通過AM1-BCC程序計算配體的原子電荷, 因此不需要提供依據RESP程序計算原子電荷所需的與高斯計算有關的所有參數. 此外, 所有配體總體上都是中性的(non_neutral_ligands=0), 因此不需要指定其中定義配體總電荷的lig_charge_file項.

運行FEW:

在上面的leap_am1文件中根據您的tutorial文件夾的位置更改了紅色的路徑之後, 就通過以下方式開始電荷計算和參數準備步驟:

perl$FEW/FEW.pl MMPBSA /home/user/tutorial/cfiles/leap_am1

在工作站機器上執行此命令可能需要大約半個小時, 因為計算原子電荷是嚴格的計算.

輸出:

當第一個設置步驟成功結束時, /home/user/tutorial目錄中會出現一個名為leap的新文件夾. leap文件夾中包含每個配體的子文件夾, 其中包括用LEaP設置MD模擬所需的結構和參數文件.

注意: 在本教程中, 演示了最簡單的過程, 即使用AM1-BCC程序確定原子電荷的模擬準備. 如果您希望使用RESP電荷設置模擬, 請參閱手冊並使用示例文件$FEW/examples/ command_files/commonMDsetup/leap_resp_step1和$FEW/examples/command_files/ commonMDsetup/leap_resp_step2.

B: MD輸入文件的準備

MD模擬的輸入文件可以使用下面展示的命令文件setup_am1_3trj_MDs來準備. 有字元的顏色與上面的步驟1A中表示的意義相同.

輸入/輸出信息:

命令文件中定義輸入和輸出數據的參數部分等同於上述準備步驟1A中所示的部分. 無論所需的功能如何, 都必須在所有的FEW命令文件中定義此部分. 否則, FEW將不知道在哪裡查找和如何處理輸入文件以及在哪裡輸出結果.

MD模擬的設置:

模擬是依據三軌跡的方法使用AM1-BCC電荷來準備. 已知存在於因子Xa蛋白質中的二硫鍵在SSbond_file項指定的disulfide_bridges.txt文件中定義, 其形式為由製表符(TAB鍵)分割的二硫鍵的殘基對的列表. additional_library文件中則提供了受體中結合的鈣離子參數. 其他參數無需設置, 因為ff12SB力場中有可用的鈣離子參數, 在受體結構中將使用力場中的這些默認值.

平衡: 使用equi文件夾中提供的模板文件準備平衡. 在此使用的平衡過程中, 分子體系首先進行兩次連續的最小化, 依次對溶質分子(即受體和/或配體)先使用強(minntr_h.in)約束再使用弱(min_ntr_l.in)約束. 在受體部分, 從殘基1到no_of_rec_residues的所有殘基被約束. 在最小化之後, 執行400 ps對溶質弱約束的MD平衡. 首先在NVT條件(md_nvt_ntr.in)下將溫度升高到300 K, 然後在NPT模擬(md_npt_ntr.in)中將密度調整到1g/cm3, 最後在50 ps NVT模擬(md_nvt_red.in)中逐漸消除對溶質的約束.

MD輸入文件的順序在MDequil_batch_template項指定的批處理腳本中定義. 如果要使用不同的平衡步驟, 則需要在批處理文件中調整平衡模板文件和sander程序的調用. 當然, 這裡使用的平衡步驟一般都應該將分子體系完全平衡. (在$FEW/examples/input_info/equi目錄下提供了本FEW教程中這個平衡過程的模板文件). 在MDequil_batch_template腳本中有關計算環境的一些具體設置通常需要根據本地計算機的進行調整. 如果您要使用上述的平衡過程, 請僅更改腳本的Fix variables聲明前的第一部分,

生產: 用於MD生產步驟的輸入文件是基於模板文件MD_prod.in創建的, 因此總生產時間為2 ns. 請注意, 這裡短的模擬時間僅供演示之用. 通常需要更長的模擬時間才能獲得完全平衡狀態的軌跡, 並從中提取具有代表性的快照. 平衡階段的最終坐標文件被用作生產階段的輸入文件, 基本文件名由restart_file_for_MDprod設置為md_nvt_red_06. 由MDprod_batch_template項指定的, 用於MD生產階段的批處理模板文件中的Fix variables聲明前和Re-queue部分需要根據本地計算機中運行MD模擬的設置調整.

運行FEW工具

當確認正確指定了基本輸入和輸出路徑後, 就可以調用FEW工具來設置MD模擬了:

perl $FEW/FEW.pl MMPBSA /home/user/tutorial/cfiles/setup_am1_3trj_MDs

輸出:

當FEW成功運行完成後, 在基本輸入/輸出目錄(即, tutorial文件夾)中會出現一個名為MD_am1(圖2)的新文件夾. 正如MD文件夾(MD-am1)的名稱已經交待的那樣, 它包含使用AM1-BCC電荷運行MD模擬的文件. 當您進入MD_am1文件夾, 您會看到每個配體的子文件夾和一個受體的子文件夾(rec). 每個子文件夾中都有一個包含初始坐標和拓撲文件的Cryst文件夾, 以及com(複合物)和lig(配體)文件夾. 後兩個文件夾中又包括各自體系中的equi(MD平衡相的所有文件)和prod(MD生產相的所有文件)兩個文件夾.

圖2: 在步驟1B中由FEW創建的MD文件夾結構示意圖

現在, 調用MD_am1文件夾中的shell腳本qsub_equi.sh來開始MD模擬.

/home/user/tutorial/MD_am1/qsub_equi.sh

當完成了所有配體/複合物的平衡階段, 就可以通過調用MD_am1文件夾中相應的shell腳本來開始生產階段.

/home/user/tutorial/MD_am1/qsub_MD.sh

注意: 如果您不想運行MD模擬, 您可以將$FEW/examples/tutorial中的整個MD_am1文件夾複製到你的tutorial文件夾下.

利用得到的MD模擬的軌跡, 可以繼續設置後面的自由能計算.

第二節: MM-PB(GB)SA 計算

在MM-PBSA和MM-GBSA方法中, FEW都支持使用隱式溶劑分子力學計算結合自由能. 設置MM-PB(GB)SA計算的先決條件是在MD模擬文件夾中已經有MD軌跡, 本例中通過FEW的MD設置功能生成這些軌跡(參見本教程的步驟1A和1B). FEW允許通過1-軌跡和3-軌跡方法準備MM-PB(GB)SA計算[1]. 對於在本教程第一節(表1)中已經準備好的三種因子Xa抑製劑的MD模擬, 這兩種方法在這裡都將作為可選方案.

有關MM-PB(GB)SA方法的基本介紹, 請參閱文獻和高級教程A3. 在教程A3中闡述的結合自由能ΔGbind, solv計算的基本原理也適用於本教程中進行的計算. 然而, 由於假定配體結合時構型熵的變化沒有顯著差異, 所以不考慮熵對結合的貢獻. 故, 我們將這裡通過FEW工具依據MM-PB(GB)SA方法計算得到的結合能稱為有效結合能, ΔGeffective.

請注意: 這裡使用的MM-PB(GB)SA計算樣本是基於本教程第一部分準備的2 ns MD模擬的快照. 因此, 不能期望用於計算的快照代表完全平衡的結構. 此外, 所考慮的快照數量較少, 進一步增加了預測的不確定性. 故, 這裡展示的MM-PB(GB)SA計算只應被看作是演示FEW功能的一個例子. 對於任何「現實生活」的研究, 強烈建議徹底分析MD軌跡, 以確定用於MM-PB(GB)SA計算的代表性結構集, 因為只有基於這樣的集合進行的MM-PB(GB)SA計算才可得到精確的結合能預測. 有關平衡狀態考察的一些基本分析可以在教程A3第2節中找到.

A: 3-軌跡方法

由於我們已經根據第一節的3-軌跡方法設置和運行了MD模擬, 所以我們可以直接從MM-GBSA的設置開始, 3-軌跡方法的計算使用FEW命令文件.

這個命令文件由以下幾個部分組成:

輸入和輸出目錄/文件的位置:

與本教程第1節中討論的輸入/輸出信息部分相同.

MM-PBSA / MM-GBSA計算設置的一般參數:

使用由AM1-BCC方法確定的配體原子電荷產生MD軌跡, 並根據3-軌跡方法準備MM-GBSA計算(1_or_3_traj=3). 為了在沒有水的情況下設置拓撲文件, 再次使用鈣離子的附加庫文件, 並且使用位於本地AMBER安裝目錄中的bin文件夾下的mm_pbsa.pl腳本進行MM-GBSA分析.

坐標(快照)提取的參數:

請求在從本教程的第1節生成的2 ns軌跡快照中沒有間隔(offset_snapshots=1)的提取1到100幀的鏡像快照.

MM-PBSA / MM-GBSA分析:

根據Onufriev[6]等人的方法計算MM-GBSA, 相應的, 在AMBER12中需要設置igb=2. 結構上結合的離子被認為是受體的一部分, 因此受體殘基的總數達290. 計算快照51至100的部分, 間隔為1.

計算參數設置為將/home/user/tutorial目錄作為基本輸入/輸出目錄以串列方式運行mm_pbsa.pl. 批處理腳本MMPBSA.sge作為作業提交的模板批處理文件. 請根據您的計算環境配置修改此腳本中Prepare calculation前面的部分.

運行FEW:

如果確保在命令文件中正確指定了基本輸入/輸出目錄的路徑, 則可以使用以下命令開始基於FEW工具MM-GBSA計算的設置:

$FEW/FEW.plMMGBSA /home/user/tutorial/cfiles/mmpbsa_am1_3trj_pb0_gb2

輸出:

在成功完成FEW運行之後, 在基本輸入/輸出目錄(即tutorial目錄)中會出現一個名為calc_a_3t的新文件夾. 該文件夾(圖3)根據所採用的方法命名(a=am1, 3t=3-trajectory approach), 其中包含每個配體的子文件夾, 配體文件夾中包含: 不含水的拓撲(文件夾topo), 提取的快照 (文件夾snapshots)以及MM-GBSA計算的輸入文件(文件夾s51_100_1).

圖3: 由FEW創建的MM-PB(GB)SA計算的文件夾結構示意圖.

執行位於calc_a_3t文件夾中的qsub_s51_100_1_pb0_gb2.sh腳本以啟動MM-GBSA計算.

計算的結果文件將會出現在文件夾中

因為如果研究十個或更多的配體, 查找所有配體的各自的結果會是一個大規模的工作, 所以FEW中提供了一個允許自動提取最終ΔGbind, effective和各自ΔE貢獻的腳本. 該腳本需要一個文本文件, 用來指定需要考慮的配體的名稱(每行一個名稱), 如:

要執行該腳本, 請進入文件夾並通過以下方式調用腳本:

最後一項由「_」分隔, 分別對應於first_PB_snapshot, last_PB_snapshot和offset_PB_snapshots. 運行結束後, 您將在/home/user/tutorial/calc_a_3t文件夾中找到一個名為pb0_gb2.txt的文件, 裡面匯總了最終結果(請參見下文).

在該文件中給出的能量貢獻對應於配體結合時的能量變化, 即複合物的形成過程, 分別為:

靜電能量

范德華相互作用

非極性溶劑自由能

極性溶劑自由能

總結合能, 即ΔGeffective

計算出的相對結合能顯示了預期的趨勢. 儘管平均有效結合能的發現是有希望的, 但是我們不應該忘記在MM-GBSA計算中考慮的構象集合的尺度之小(僅考慮來自2ns MD模擬的50個快照). 當你在文件/home/user/tutorial/calc_a_3t/L51c/s51_100_1/pb0_gb2/L51c_statistics.out.snap中的DELTA部分尋找例子時, 您會發現為單個快照計算的ΔGeffective值可以從約-165.6至+84.07 kcal/mol. 因此標準偏差很高(60.8 kcal/mol). 計算能量大幅度波動的一個重要因素是內能中的噪音. 在3-軌跡方法中, 複合物, 受體和配體的內能不會相互抵消, 因此即使結構之間的微小差異也可以對計算出的結合能量產生很大的影響. 在1-軌跡方法中避免了這種噪音, 該方法僅根據複合物的結構計算結合能量, 不考慮配體和受體的遊離形式.

B: 1-軌跡方法

儘管事實上在所謂的1-軌跡方法中只考慮了來自複合物模擬的結構, 即未考慮受體和配體在未結合和結合前後的結構差異, 因此這種方法經常被應用並且表明具有良好的相對結合能預測能力. 這種方法的一大優點是對計算的要求較低, 因為只考慮了複合物的MD模擬. 此外, 該方法還可以消除內能的貢獻, 從而具有更低的ΔE噪音值.

因此, 由於經常使用1-軌跡方法來估計相對配體結合能, 這裡也將演示基於這種方法的MM-PB(GB)SA計算的設置. 一般步驟與上述的3-軌跡方法所示的步驟一樣. 我們可以按照步驟2A進行. 首先為FEW創建一個名為mmpbsa_am1_1trj_pb3_gb0的命令文件, 它與上面展示的mmpbsa_am1_3trj_pb0_gb2文件一致. 然後設置1_or_3_traj=1, GB=0和PB=3. 因為這裡將執行在計算上更費時的泊松-玻爾茲曼計算, 所以可通過設置first_PB_snapshot=81來減少考慮的快照的數量以節省計算時間. 如果您現在使用此命令文件運行FEW, 則會根據1-軌跡方法獲取帶有解析半徑的MM-PBSA計算的輸入文件. 計算的最終結果如下所示. 儘管模擬時間短, 所考慮的快照數量較少, 並且通過1-軌跡方法進行的近似, 但計算的結合能ΔGeffective(E_TOT)依然顯示了這些配體的預期趨勢(參見表1).

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！