以熱力學積分分子動力學模擬方法計算細胞外信號調節激酶與其抑制劑的結合自由能

全文

(1)國立臺灣師範大學化學系碩士論文 Department (Graduate institute) of Chemistry,. National Taiwan Normal University Master Thesis. 指導教授：孫英傑博士 Advisor: Ying-Chieh Sun, Ph.D.. 以熱力學積分分子動力學模擬方法計算細胞外信號調節激酶與其抑制劑的結合自由能. Computation of Binding Free Energy of Inhibitor-Erk Kinase Complexes Using Thermodynamic Integration Molecular Dynamics Simulation. 研究生:陳博晉 Graduate: Po-Chin Chen. 中華民國 101 年 6 月 June, 2012.

(2) 謝誌時光飛逝，歲月如梭，光陰似箭，時光荏苒。碩士班這兩年的時間很快的過去了。遙想當年大學畢業時青澀模樣，已成為記憶中的一部分了。如今的自己相較於過去成長了許多，未來的挑戰也正等著我，我將義無反顧、勇往直前、全力以赴。回想這兩年來在師大生活的點點滴滴於畢業過後都將成為過去，學生的求學階段永遠是最好的，能夠在一個自由的校園環境下自主學習成長以及向同學請教、學習，這於人生長河中都是非常寶貴且值得珍藏回憶的過程。除此之外在校修課、實驗室做研究、出外參與年會都是一段美妙的經歷。這兩年當中很感謝我的指導教授孫英傑博士的諄諄教誨，不管是在學業上或者是待人處事方面老師都給了我許多具體建議，使我能將自己的缺點加以改進，讓自己比以前更好。另外也感謝能李冠群教授以及楊小青教授撥冗擔任學生的口試委員，給予研究上的建議與指導。此外我也很感謝我師大的同學(不論是以前的同學或是進入以後才認識的同學)，有了你們的支持鼓勵以及分享使我在這兩年的碩士班生涯遇到困難或挫折能夠坦然面對、虛心接受、妥當處理、安心放下。實驗室當中有一群好相處以及熱誠的夥伴們也使我們實驗室充滿歡笑，學長昭同、書吉、宗儒，在我剛來實驗室有很多知識不懂之時，有了你們的指導使我熟悉這個環境且學會專業上的知識。同學冠緯、孟璇、永富、彥樓也謝謝你們陪我度過這兩年，不論是課業或是日常生活中得到你們的幫忙以及分享，我的心中充滿感激。除此之外，興洲、育劭還有嘉仁，也謝謝你們跟我討論分享使我獲益良多。還有實驗室來自中國大陸的博士後研究員王軍學長，給予我許多做科學研究上的指導以及建議。在人生的旅途中因為有你們讓我更加難以忘懷在師大的生活。很高興認識你們，天下沒有不散的筵席，希望以後能有緣和你們敘舊。另外感謝國家網路高速計算中心所提供的計算資源；本研究若沒有國網中心的支持，研究的內容必將無法完整。當我在國網中心計算碰到問題請教何博士，何博士總是以最快的時間內協助我們克服困難，才得以延續本次的研究，謝謝您，何博士。最後也要很感謝我的家人，我的父母親對於我想要做的事情都採取鼓勵的態度，支持孩子完成自己的夢想。生活中給予:高關懷、多監督、少責備的態度，因此我有足夠的空間追逐自己的理想完成自己的抱負，我心存感恩。最後我也感謝我妹，讓我在家除了爸媽之外還有多一個人可以互動。師大之旅，畢生難忘。 i.

(3) 中文摘要 Erk 激酶是 Ras/Raf/MEK/ERK 訊號傳遞路徑的下游蛋白質，此路徑在某些癌症當中失調。我們選擇此蛋白質激酶並且使用分子動力學模擬來設計此蛋白質激酶新的抑制劑。先前的計算當中發現熱力學積分分子動力學模擬方法於計算蛋白質激酶一組抑制劑得到的相對自由能與實驗結果相近。在目前的研究上我們對新的類似化合物進行 TI-MD 計算是為了要搜尋出更好的抑制劑。一個在苯環上帶有 OH 官能基的類似物被發現比 82A 有更好的親和性約 1.1 千卡/莫耳。除此之外，我們探討於 TI-MD 方法中一步變換與三步變換上的表現之比較。另外，也討論關於不同的 λ 數目在計算上對於相對自由能上的影響。這些結果應該會對於未來 TI-MD 模擬以及這種類的抑制劑之設計會有一些幫助。. 關鍵字:分子動力學模擬、自由能、細胞外信號調節激酶、熱力學積分。. ii.

(4) ABSTRACT Erk is the downstream protein in the Ras/Raf/MEK/ERK signal transduction pathway which is deregulated in many cancers. In the present study, we select this kinase to design new inhibitors for this kinase using MD simulation. Previous computations have shown that thermodynamic integration (TI) MD simulation method is able to give relative binding free energy for a pair of inhibitors for this kinase in good agreement with experimental results. In the present study, we carried out TI-MD computations for new analogous compounds in searching for better inhibitors. An analog with –OH functional group on the benzene ring was found to have better affinity than the 82A ligand by 1.1 kcal/mol. In addition, we investigate how single mutational step performs compared with the triple mutational steps in the TI-MD method. Furthermore, effect of number of lambda points in the computed relative binding free energy was investigated as well. These results should be useful for further TI-MD simulation and inhibitor design of this kind.. Key word : Molecular Dynamics Simulation、Free Energy 、 Extracellular signal-regulated kinases(ERK)、Thermodynamic Integration(TI) iii.

(5) 目錄(CONTENTS) 口試委員會審定書 ........................................................................................................... # 謝誌. i. 中文摘要. ii. ABSTRACT iii 目錄(CONTENTS) ...........................................................................................................iv 圖表目錄 (LIST OF FIGURES) .....................................................................................vi 表格目錄(LIST OF TABLES) ...................................................................................... viii Chapter 1. 緒論 (Introduction) ................................................................................... 1. 1.1. 前言 (Preface) .......................................................................................... 1 1.1.1 藥物開發(Drug development) ........................................................................ 1 ERK 蛋白質激酶與癌症的關聯 (The Relationship Between ERK. 1.2. Protein Kinase and Cancer) .................................................................... 4 1.3. 研究目標(Research Objectives) .............................................................. 8. Chapter 2. 實驗方法 (Experimental Methods) ........................................................ 10. 2.1. 實驗理論 (Theories) .............................................................................. 10 2.1.1 分子動力學模擬 (Molecular Dynamics Simulation) ................................... 10 2.1.2 分子動力學模擬理論 (Molecular Dynamics simulation theory) ........................................................................................................... 10 2.1.3 參數組 (Parameter set) ................................................................................. 12 2.1.4 應用程式 (Application) ................................................................................ 14 2.1.5 熱力學積分 (Thermodynamic Integration, TI) ............................................ 15. 2.2. 實驗設定 (Modeling Settings) .............................................................. 18 2.2.1 分子動力學模擬設定 (Molecular Dynamics Simulation settings ) ........................................................................................................ 18 2.2.2 RESP (Restrained electrostatic potential) ..................................................... 19 2.2.3 Functional group transformation .................................................................. 20 2.2.4 Cheng-Prusoff equation ................................................................................ 20 2.2.5 TTA v.s. STA................................................................................................. 23 iv.

(6) 2.2.6 Comparison of different λ points .................................................................. 24 Chapter 3. 結果與討論 (Results and Discussions) .................................................. 25. 3.1. TTA (Three transformation approach) relative binding free energy( λ points = 9) ............................................................................................... 25 3.1.1 TTA：CFO to 82A (solution phase) ............................................................. 25 3.1.2 TTA：CFO to 82A (complex phase) ............................................................ 27 3.1.3 Summary of CFO to 82A .............................................................................. 32. 3.2. Comparison of TTA with STA .............................................................. 33 3.2.1 Results of STA (CFO to 82A) ....................................................................... 33 Result of different λ points in TTA ....................................................... 39. 3.3. 3.3.1 λ points equal to 9 ......................................................................................... 39 3.3.2 λ points equal to 5 ......................................................................................... 42 3.3.3 λ points equal to 4 ......................................................................................... 45 3.4. Difference between each lambda points ............................................... 49. Chapter 4. 結論 (Concluding Remark)..................................................................... 52. 英文專有名詞中文翻譯 ................................................................................................. 54 參考文獻(Reference) ...................................................................................................... 56. v.

(7) 圖表目錄 (LIST OF FIGURES) 3. 計算精準度與時間的比較 ...................................................................................................... 2. 11. Raf/MEK/Erk 路徑概觀 ......................................................................................................... 4. 圖表 1-1 圖表 1-2. 圖表 1-3 信號傳遞路徑 .......................................................................................................................... 5 圖表 1-4. ERK2 以緞帶顯示之結構 .......................................................................................................... 6. 19. ERK2 結構分析示意圖 ........................................................................................................... 6. 20. 圖表 1-6. 尚未磷酸化之 ERK2 結構 ...................................................................................................... 7. 圖表 1-7. 82A 結構圖............................................................................................................................... 8. 圖表 1-5. 25. 原子型態表.......................................................................................................................... 13. 25. MD 模擬的基本流程圖表 .................................................................................................... 14. 26. 熱力學循環計算不同配體下的相對自由能 ........................................................................ 15. 圖表 2-4. 10a. 各種類似化合物之活性 ..................................................................................................... 22. 圖表 3-1. λ=0 CFO 進行變換形成 λ=1 82A ............................................................................................. 25. 圖表 2-1 圖表 2-2 圖表 2-3. 圖表 3-2 配體 λ=0 未含水結構:CFO. 圖表 3-3 配體 λ=0 含水結構:CFO ........................................ 25. 圖表 3-4 配體 λ=1 未含水結構:82A. 圖表 3-5 配體 λ=1 之含水結構:82A..................................... 26. 圖表 3-6. 2OJJ 緞帶形狀結構紅色部分為 Helix 部分紫色為 strand 部分 ......................................... 27. 圖表 3-7. 2OJJ 加入水盒子後之構型。 ................................................................................................. 28. 圖表 3-8. CFO 以及 82A 在 complex phase 下重疊示意圖 .................................................................... 28. 圖表 3-9. 82A 與周圍胺基酸交互作用情形 .......................................................................................... 31. 圖表 3-10. CFO 與周圍環境胺基酸交互作用情形 ................................................................................. 31. 圖表 3-11. solution phase 9 個積分點與 19 個積分點 .......................................................................... 34. 圖表 3-12. STA complex phase lambda points=9 .................................................................................... 35. 圖表 3-13. solution phase 下 TTA 與 STA 及調整參數和不同 λ 數目比較 ............................................. 35. 圖表 3-14. TI2 complex phase(0-0.6ns) λ points=9 ................................................................................. 39. 圖表 3-15. TI2 complex phases(o.6-1.2ns) λ points=9 ............................................................................ 40. 圖表 3-16. complex phase(0.6-1.2ns) λ points=9 .................................................................................... 40. 圖表 3-17. solution phase (0-0.6ns) λ points=9 ...................................................................................... 41. 圖表 3-18. TI2 complex phase(0-0.6ns) λ points=5 ................................................................................. 42. 圖表 3-19. TI2 complex phase(0.6-1.2ns) λ points=5 .............................................................................. 43. 圖表 3-20. complex phase(0.6-1.2ns) λ points=5 .................................................................................... 43. 圖表 3-21. solution phase (0-0.6ns) λ points=5 ...................................................................................... 44. 圖表 3-22. complex phase(0.6-1.2ns) λ points=4 .................................................................................... 46. 圖表 3-23. TI2 complex phase(0-0.6ns) λ points=4 ................................................................................. 46 vi.

(8) 圖表 3-24. TI2complex phase (0.6ns-1.2ns) λ points=4 .......................................................................... 47. 圖表 3-25. solution phase (0-0.6ns) λ points=4 ...................................................................................... 48. vii.

(9) 表格目錄(LIST OF TABLES) 表格 1-1. 82A 於對不同蛋白質之 Ki 值 ................................................................................................... 8. 表格 3-1. solution phase Total Gibbs energy .......................................................................................... 26. 表格 3-2 complex phase TI1&TI3 free energy ........................................................................................ 29 表格 3-3 complex phase TI 2 使用 6 個亂數平均而得之能量 .............................................................. 29 表格 3-4. TI solution & complex phase of total free energy.................................................................... 32. 表格 3-5. STA(production run 600ps)計算結果 ...................................................................................... 33. 表格 3-6 更改 scbeta 數值為 7 後之比較 ............................................................................................. 36 表格 3-7. TTA(solution phase)電荷減半所得之自由能 .......................................................................... 37. 表格 3-8. STA(solution phase)電荷減半所得之自由能 .......................................................................... 37. 表格 3-9. TTA(solution phase)電荷減半與未減半之比較 ...................................................................... 38. 表格 3-10. STA(solution phase)電荷減半與未減半之比較 .................................................................... 38. 表格 3-11. λ points=5 能量總整理 ........................................................................................................ 42. 表格 3-12. λ points=4 能量總整理 ........................................................................................................ 45. 表格 3-13. solution phase at different lambda points ............................................................................ 49. 表格 3-14. complex phase at different lambda points ............................................................................ 50. 表格 3-15 不同 lambda 點所計算出之相對自由能差值...................................................................... 51. viii.

(10) Chapter 1. 緒論 (Introduction). 前言 (Preface). 1.1. 細胞訊息傳遞與代謝是人類維持生理運作的重要功能，其中一部份訊息傳遞過程是由細胞膜上的蛋白質激酶做為訊息接收者，接收膜上的蛋白質激酶做為訊息接收者接收膜外的訊息，蛋白質激酶活化後把訊息從膜外傳遞到細胞內，再由細胞質內特定蛋白質激酶經由磷酸化方式把訊息傳遞到下游，最後將訊息傳遞到核中決定核內的 DNA 基因表現來調控細胞功能，例如:細胞的生長、轉移以及自我凋零。在化學或者是生物化學領域中，蛋白質激酶為一種酵素，它又稱做為 phosphotransferase，在生物體當中三磷酸腺苷會和蛋白質激酶結合以進行磷酸化，調控下游蛋白質的活性，而三磷酸腺苷為一種能量分子，在與蛋白質激酶進行磷酸化作用後，三磷酸腺苷會轉移一個磷酸根到受體蛋白質。有些時候若某一個蛋白質失調過度表達或者突變有可能會導致細胞異常的增生或者是變異而產生癌細胞，起因為傳遞路徑上一個或多個蛋白質變異所導致，因此當激酶是標靶時，我們藥物設計的策略就著重在找出一個小分子化合物來和三磷酸腺苷競爭，搶在蛋白質激酶和三磷酸腺苷之前結合，就可以使得細胞傳遞訊號的途徑被抑制，來達到降低活性的效果。也就是說需找出一個和蛋白質結合能力更強的小分子，以達到治療癌症的效果。. 1.1.1. 藥物開發(Drug development). 目前在開發藥物的成本依舊居高不下，投入大約 500 億美金僅能成功開發出約 20 顆新藥 1，而當務之急就在於如何以較低的成本、更快的速度開發出一顆新藥，到目前為止這依舊是一大挑戰，加以現今人口老化，若能克服以上問題，對於人類來說會是一大福音。電腦模擬輔助藥物設計在這樣的環境下對於藥物開發會是一項很重要的利器。因為它能夠節省做實驗所需購置實驗器材的成本、以及材料的成本。若能使用電腦模擬精準的得到新開發藥物對於癌症的抑制效果那會是最理想的狀態。現今新型藥物的開發，已逐漸從藥物發現轉為藥物設計。現今電腦模擬的技術還無法完全的精準預測新開發藥物抑制效果，某種程度上還是有誤差，加以存 1.

(11) 在的化合物非常之多，要用甚麼樣的方法能夠達到快速、準確的得到有潛力的藥物呢？這就是現今研究的課題。使用分子嵌合大量篩選的方法往往只能夠定性，而無法定量。若要得到較為精確定量則需使用分子動力學模擬，而計算結合能的方法有數種，本篇所採取的方法為熱力學積分得到相對自由能。癌症藥物的開發往往於一項藥品開發出來之後，須經過重重實驗和審核，才能上市。首先必須通過動物實驗，接著經過三期的人體實驗，每期人體試驗分別要測試不同的人數，所需要投資高昂的成本，成本回收期也相當的長，2 在人體實驗中遇到失敗，先前投資的金錢就無法回收，可見生醫產業的進入門檻和條件限制相當高，也因此可見如果能夠運用電腦模擬的結果來提供有潛力的化合物以進行下一步的合成實驗，最後再進入生物實驗，此將可提高藥物發展成功的可能性。故發展電腦模擬以提供藥物設計上的建議是有價值的，它能夠快速搜尋有潛力的藥物，使得研究藥物到上市的時間縮短，另外也能降低成本，提升成功機率。. 圖表 1-13. 計算精準度與時間的比較. 目前電腦模擬無法做到百分之百精準，還是會有誤差的存在。而如今存在於世界上的化學結構非常的多，我們要如何知道哪一個化合物是有潛力的藥物呢? 4 所以一般來說可以先使用大量、快速篩選的方式，找出有可能成為藥物的化合物，此即定性的方式；接著再使用較為精準的方式，較為精確的比較出個別藥物結合蛋白質強度。從上圖可以看到從下至上模擬時間以及準確性的比較：圖中最下面的 score 即為快速篩選的函數，目前世界上有多家公司在發展此類軟體，每家公司發展出的軟體準確性不盡相同，於市面上販售的價格也不相同。此端看研究人員 2.

(12) 的需求而定。它的優點為操作難度低，較容易讓使用者能夠上手，並且能夠在短時間內搜尋出許多不同種類官能基的化合物，準確定雖然不高，但也由於篩選出的化合物結構可能差異性大，對比於現今存在可能被專利過的化合物會較為不同，而能避開專利的限制。但也有其缺點：準確度較為低，只能得到大概的結果。在圖表中間以及上面為計算自由能的方法 5，以本實驗來說：計算此自由能所使用的軟體為 AMBER。用上圖的計算方法所得到的計算結果，以圖表來說越上方的方法會越為精準 6。以 MM-PBSA7 的方法來說，7b 誤差大約在 3~5 kcal/mol 左右， TI 方法誤差能夠降低到 1 kcal/mol 左右。不過此類越上方的計算方法也有其缺點以及困難，所需計算的時間相當漫長，且要有大量的高速計算主機，另外也需要花上不少的時間使得計算能量得到收斂的數值 8。加以操作上的難度較大量篩選的方法來的高，9 並不是這麼容易上手 8-10，所以研究上還是有科學家用 MM-PBSA10b，或者是 Docking 來進行化合物篩選。. 3.

(13) 1.2. ERK 蛋白質激酶與癌症的關聯 (The Relationship Between ERK Protein Kinase and Cancer). 圖表 1-211 Raf/MEK/Erk 路徑概觀上圖為 MAPK pathway 訊息傳遞路徑示意圖，細胞質上 EGF 接到受體 EGFR 活化了細胞質區域上之酪胺酸激酶活性，而後的磷酸化進行調控下游的活性，最後將訊息傳入細胞核內決定 DNA 基因表現。 ERK 的全名為：(Extracellular signal-regulated kinases 細胞外信號調節激酶)。 12. ERK 在細胞增殖中之信號傳遞扮演了重要的角色。研究 ERK 是研究哺乳類. MAPK 路徑當中最好的研究。約三分之一癌症被發現與 ERK 路徑失調有關，此一途徑控制了細胞基本進程：包括細胞增生、存活、分化、運動。13 以歷史上來看， ERK 即代表著細胞增殖，現階段的研究看到 ERK 的失調與某些腫瘤型態息息相關；包括：結腸癌、腎癌、胰腺癌、卵巢癌。以下圖來看：14ERK 被 MEK 磷酸化，而 MEK 又被 RAF 磷酸化，以下為信號傳遞的順序:. 4.

(14) 圖表 1-3. 信號傳遞路徑. ERK 是由各式各樣的細胞外訊號所活化，大多數和癌症有關的病變，導致構型活化 ERK 訊號發生在這個途徑較前面的步驟，即受體酪氨酸激酶過度表達、活化受體酪氨酸激酶突變、保持活化配體分泌或旁分泌、RAS 突變或 RAF 突變。然而也有核轉錄因子失調或是放大的標靶，最顯著的為 myc 和 AP-115。 ERK 可分為 ERK1 以及 ERK2，ERK2 又被稱為 MAPK1(Mitogen-activated protein kinase 1)不論是在透過增加膜結合生長因子受體突變或者是上游蛋白質表達造成 ERK 蛋白過度表達被認為與人類多項的癌症病變有相關性。若可針對 ERK 信號來加以抑制會是有潛力的策略來防止癌細胞的擴散 16。大多數的 kinase 抑制劑缺乏特異性，因為他們是針對蛋白質激酶家族中的 ATP 結合區域。17 現今 ERK 的抑制劑發展陸陸續續在文獻當中可以搜尋找到，其中美國 vertex 公司在期刊中發表了幾篇關於 ERK 的論文，本篇論文所做的研究即依據 vertex 公司在 Journal of Medicinal Chemistry(JMC)於 2007 年所發表的文章為基礎，進一步的尋找更為理想的 ERK 抑制劑。而 vertex 公司也於 2009 年再度發表一篇關於 ERK 抑制劑的論文。18. 5.

(15) 圖表 1-4. ERK2 以緞帶顯示之結構. ERK2 之未磷酸化 3 度空間緞帶形狀之結構。此圖顯示了 ERK2 磷酸化位置以及 docking 區域之空間關係，紅色部分為 Thr183、Tyr185，他們是行磷酸化的殘基(residue)；綠色部分為 Thr157、Thr158，藍色部分為 Asp316、Asp319。藍色以及綠色的部分被稱為 docking residues。. 圖表 1-519. ERK2 結構分析示意圖. N 端為殘基 1-109 以及 320 到 358 是由β-strands(圖上以綠色表示) 、兩個 helices(C 藍色以及 aL16 洋紅色)組成。 C 端為殘基 110–319 大部分為 helices(藍色)、洋紅色 Lip 磷酸化位置。. 6.

(16) 圖表 1-620. 尚未磷酸化之 ERK2 結構. 以圖表 1-320 來說 ATP 與 ERK2 結構的位置在結合區域接面較靠近 active site(活化位)內部而蛋白質基質是接在表面，控制激酶活性的兩個殘基 Tyr185、 Thr183 位於活化迴圈(activation loop) 21，ERK2 單一殘基磷酸化並不會造成活性的增加，也不會利用酸性氨基酸來替換磷酸化位置。這可能因為是在進行磷酸化時，構性一定會跟著改變的特性。硫酸根離子所在的位置是在磷酸活性結構上的位置。而 D316、D319(aspartic acid residue)為：D 區域被建議的結合位置。. 7.

(17) 1.3. 研究目標(Research Objectives) Erk 在許多癌症的病變當中有一個共同的特色就是與 Ras/Raf/MEK/ERK 路徑. 有關聯性。而我們希望在既有的文獻當中參照並且設計對於比現有藥物分子更為理想的藥物分子。我們針對文獻上對於 ERK2 有不錯抑制效果的化合物 82A 進行了小幅度的官能基變換。並且利用 TI 的方法做計算，最後去和原先文獻上的化合物做抑制效果的比較。. 本篇研究希望能夠找出相較於 vertex 公司於 2007 年發表論文所做的一系列化合物中能夠找到更好的抑制劑 10a。其中那篇論文當中有報告一解出的結晶結構(pdb code: 2OJJ)。從蛋白質資料庫下載的結構 2OJJ 所得到的資料顯示：結晶結構當中的配體為：82A. 圖表 1-7. 82A 結構圖. 從那篇 paper 當中可發現 82A 對於 ERK2 有相當好的抑制效果如以下所示: 表格 1-1. 82A. 82A 於對不同蛋白質之 Ki 值 ERK2,Ki,μm. JNK3,Ki,μm. 0.002. >4. 8.

(18) 從上表可看出對於 ERK2 而言：Ki 值可以到達 0.002μm，也就等同於 2nm，以一個新設計出來的化合物而言，單就 Ki 值來說 82A 是個還不錯的抑制劑。本篇即希望以 82A 為主體架構去進行官能基的變換，而後能夠利用 TI 熱力學積分的方法計算得到經由變換後得到之新化合物相對自由能 10a。經過閱讀 vertex10a 公司論文後，本次研究上挑選出欲變換的官能基為 F 以及 CL，考量到周圍的環境決定變換成其他的小官能基，同時挑選較具極性的官能基，以及非極性的官能基來做替換以進行相對自由能的計算。所以此次決定將 F 置換成 CF3、CL 置換成 OH。本篇也探討 TI 計算運用不同的步數(step)所得到的差異，現今文獻上對此一課題探討的較為少，本論文當中也有做一部分這方面的探討。若能使用較少的步數對於計算上來說會較為方便，因為這可以減少計算的程序。相對自由能計算所選用的方法為 TI，依照 TI 計算上程序的不同可分為兩種：第一種程序為 TTA，另外一種程序為 STA22。 TTA 即計算採用三個步驟；STA 即計算採用一個步驟。目前相關研究上以 TTA 程序較為成熟，STA 探討的人較為少。而 TI 計算點的挑選上當然越多越好，這可以使得整個變換過程的能量計算相對來說得到較為精確的結果，但是每多挑選一個計算的點就會增加許多的計算量；即增加電腦 service unit(SU)的使用量，亦即增加計算的時間。所以必須要取得計算精準以及計算時間上的平衡，故本次研究參考了瑞典研究團隊挑選點上的策略 23 ，以進行不同點上所得之數值比較。另外也針對 STA 計算出來的數值與 TTA 的數值進行比較，目前探討兩者計算差異的論文不多，所以本篇也有去探討 STA 與 TTA 自由能以及相對自由能的差異，以及原因。最後也對於不同的 λ 點去進行比較，一般來說如果用比較多的 λ 點，積分出來的數值會更為接近實際的真實數值，但是用 TI 方法計算量很大，需要用到很多的計算資源，故能夠用比較少的 λ 點就能得到接近的相對自由能數值會更佳，因此希望能找到數值精準以及減少計算之間的平衡點，以期能夠達到效益最佳化。. 9.

(19) Chapter 2. 實驗方法 (Experimental Methods). 實驗理論 (Theories). 2.1. 2.1.1 分子動力學模擬 (Molecular Dynamics Simulation) 本實驗使用 Amber11 程式集計算分子動力學模擬以及小分子與蛋白質激酶間的作用力。Amber 直接翻譯成中文之意為:琥珀。但是在分子動力學模擬的領域上， Amber 是一個著名的軟體，它的全名為:Assisted Model Building with Energy Refinement。它分為兩部分，第一部分為:模擬生化分子的分子力場。另外一部分為:分子動力學模擬軟體的組件以模擬力場。起初 Amber 是由 Peter Kollman(美國加州大學舊金山分校教授)所帶領的研究團隊發展出之軟體，而現今有一個專門維護此軟體的團隊做更新以及解決程式上尚有缺陷的部分。. 2.1.2 分子動力學模擬理論 (Molecular Dynamics simulation theory) 現代理論計算化學中，分子動力學模擬是探討系統平衡與動力學性質的重要方法，利用古典分子動力學模擬為基礎，了解分子系統中原子的微觀運動，以數值方法解牛頓方程式:. F= m = 計算分子間能量之時，會建立一個力場描述分子中:鍵的伸縮、鍵角、扭角、凡得瓦作用力以及靜電作用力。24. V= +. -(. }+. 前三項所表示為:鍵結作用力；其中 k 為常數，r 為原子間距離，為角度、. 為扭. 動能障、n 為週期、為扭角、為相位。後兩項為:非鍵結作用力；當中為凡得瓦位能井深度、. 為凡得瓦半徑、代表圓周率。. 分子模擬過程中，描述分子系統的運動軌跡，是將每個原子視為一個質點，質點上有 x、y、z 三個方向，分別去解 x、y、z 方向的牛頓方程式。解牛頓方程式 10.

(20) 中需要原子的速度以及初始位置，模擬溫度當中會程式會利用亂數產生波茲曼分布來得到起始速度，初始位置、位能函數，便可以解牛頓方程式。為了得到下一個時間點，分子運動過後的位置，必須要給予時間梯步△t，經由計算後所得下一步位置以及速度，在重複此過程 N 次後的△t 後，可以得到我們所要模擬分子運動時間的總軌跡，以及分子在系統中的運動過程。. 分子動力學模擬過程當中為了模擬蛋白質真實常溫常壓下的環境，模擬計算會有許多方法去控制溫度壓力。在控制壓力項時，以改變體積大小去平衡瞬時壓力(. )。. =. a=-. 控制溫度項部分，N 個粒子的系統中，瞬時動能(. )會因為每一個△t 瞬時. 速度不同，瞬時速度不同會造成瞬間動能以及溫度改變，控制某一溫度時，會採取改變速度的方式達成，其中. 為瞬時溫度。. =. =. 但因電腦的計算時間有限，必須在短時間內得到所需要的平衡溫度，且對每個粒子達到平均升溫的效果，使系統快速且有效的達到平衡，在模擬的過程當中會加入 Langevin dynamic 方法使之平衡，此方式在改變瞬時動能時會加入 random force 以及摩擦係數兩項要素作為修正。. =m. 為摩擦係數，單位為 p. ，. -. (t). 為摩擦力， (t)為 random force，此項須. 符合高斯分布。生化分子的模擬系統中，在解牛頓方程式時，1fs 的△t 可獲得準確的運動軌跡，但是對於高頻率的運動牛頓方程式解出來的數值誤差相較於低頻率的運動為 11.

(21) 大，因此利用外加限制去控制分子運動中的某些維度，SHAKE 是一種 constraint algorithm 的方法，用來確保主體加上限制(constraint)後，在牛頓方程式下依舊能進行運算，於此利用固定重原子接上氫原子上的鍵長，此鍵結仍會有高頻率的運動，當將此鍵結加以限制，△t 由 1fs 至 2fs 時，模擬仍然可以得到準確的運動軌跡。另外在分子動力學模擬當中，邊界條件是用以外顯性水模型時，是以原立方體水盒子於各邊界複製相同水合子，如此模擬可以消除邊界的效應，也對於電腦記憶體不會造成過大的負擔。另外在生化模擬當中當系統帶有離子作用力時，為了讓電荷作用能量得到良好的收斂，於模擬時通常保持系統為電中性，簡單的方式為加入 Na+離子或者是 Cl-離子中和，Na+或者是 Cl-離子在此稱為 counterions。. 2.1.3 參數組 (Parameter set) 我們使用力場的時候，需要針對不同的模擬力場給予參數，例如:力常數、平衡的鍵長、鍵角、電荷…等等。在 Amber 軟體中有許多的參數組，每一個參數組都有各自的特色，以及名字，它會提供特定的參數給不同類型的原子。以下為 Amber 力場中的其中三組： 1. ff99 :這組提供胜肽、蛋白質、核酸的相關參數。名字的開始為 ff，後面兩個數字為年份。 2. GAFF(General AMBER force field):此組提供小分子相關參數，以促進小分子藥物在蛋白質當中的模擬，且能夠建立起非基本胺基酸的力場。 3. GYLCAM: GYLCAM 力場是用來模擬碳水化合物。. 下頁的圖表即代表 Amber 系統中對於各式樣鍵結的原子描述:. 12.

(22) 圖表 2-125. 原子型態表 13.

(23) 2.1.4 應用程式 (Application) AMBER 程式集當中有許多應用程式，以下僅介紹幾種常用的應用程式 1.Leap:用來準備模擬程式的輸入檔案 2.Antechamber:若系統中含有非標準胺基酸或者是蛋白質時可使用 Antechamber 利用 GAFF 自動參數化小分子。 3.SANDER(Simulated Annealing with NMR-Derived Energy Restraints) :為基本能量極小化和分子動力學模擬的應用程式。程式也會藉由將原子移往能量梯度較低的地方直到較低的平均梯度以達成鬆弛結構。 4.pmemd(particle Mesh Ewald Molecular Dynamics)為 sander 的版本之一，計算速度較 sander 來的快特別是在超過 8-16 個處理器進行運算時效果更為明顯 5.ptraj:提供進行 MD 模擬所產生的軌跡檔案做數值分析，例如:可分析氫鍵、鍵長、鍵角…等等。 6.MM-PBSA:是一個 script 可用來擷取 MD 模擬當中的軌跡資訊，以此資訊的結構來進行能量計算. 圖表 2-225. MD 模擬的基本流程圖表. 起始的結構可以用 LEaP 補氫補電荷加入水分子，antechamber 能夠轉檔處理一些關於 ligand 的相關資訊，AMBER 程式集本身無法判別非胺基酸的訊息，我們必須要使用到 antechamber 做為輔助建立其 force field 使得 AMBER 系統在建立起 14.

(24) 使模擬結構 prmtop 以及 inpcrd 檔案能夠包含有 ligand 的訊息。接著使用 sander 或者是 pmemd 來進行模擬，最後產生的檔案，可以利用 ptraj 來分析整段過程的…ex:RMSD，H-bond…；或是用 mm-pbsa 計算出結合能量。. 2.1.5 熱力學積分 (Thermodynamic Integration, TI). 圖表 2-326. 熱力學循環計算不同配體下的相對自由能. 一般說來直接計算蛋白質與配體之結合能要得到精確的數值較為困難，也因此我們使用一個間接的方式：去計算配體在水溶液中與在蛋白質間的自由能差。在一個已知配體 Ki 值之條件下，可以用熱力學循環去計算並且回推我們所設計化合物之 Ki 值。熱力學積分可用來計算相對自由能，以上圖為例:P 代表蛋白質，X 以及 Y 分別代表不同的配體。熱力學積分計算出數值的公式為:. ΔΔ. =Δ. Δ. = Δ. -Δ. 我們可以利用計算在溶液中配體變換的自由能與蛋白質配體的自由能得到相對自由能差，來推得蛋白質與配體 X 與蛋白質與配體 Y 之相對自由能差。如以下公式: 15.

(25) Δ. 其中. =. 代表在之下的位能函數，計算熱力學積分之時，我們必須要計算從. =0 到 =1 之間的位能，而後再以積分的方式計算出位能曲線下的面積，這面積就是我們最後要得到之相對自由能ΔΔ. 在熱力學積分中計算位能函數. 。. 時，我們使用到 softcore potentials，softcore. potentials 是 TI 計算中其中的一種 potentials，此系統的建立被簡化成 =0 與 =1。 =0 代表原子消失 =1 代表新的原子出現。使用 softcore potential 時，原子消失與新原子生成同時進行，. 代表初始的位能函數；而. 代表最終的位能函數。而其中 =0. 到 =1 歷經了配體的轉換，在準備 =0 與 =1 之 topology 檔案以及 coordinate 檔案時需使起始結構與最終結構兩者相同的原子之位置相同，而整體的結構也須合理。除了相同的原子之外，其他的特定原子在此將之視為 soft core atoms。使用 softcore potentials27 進行 transform 可避免掉 singularity effect，因此可使用 linear mixing 來計算出最後的 V(λ)。. 模擬當中為了要計算出最後的相對自由能，我們會將不同 λ 之下所算出的各個位能下面積去做積分，要算出 V(λ)我們使用到以下的 mixed potential function:. 16.

(26) 為起始哈密頓算，為擾動哈密頓算符。此程式也計算 < > ，在以上公式中， < >= 的數值就是從. ，同時這也被稱做為 linear mixing。最後我們積分出九個點所得到到. 這段期間的自由能。. 17.

(27) 2.2. 實驗設定 (Modeling Settings). 2.2.1. 分子動力學模擬設定 (Molecular Dynamics Simulation settings ). 模擬當中所要使用的起始結晶結構是由蛋白質資料庫(RSCB protein data bank) 所下載而得，我們的模擬狀態都是在水溶液的條件下，在模擬前我們會用 LEaP 創造水溶液的環境，過程為先將蛋白質結晶結構當中的結晶水取出而後加入一個立方體的水盒子(TIP3P)。整段模擬的過程分為以下四段:1.建立起使結構 2.能量極小化 3.升溫 4.等體積模擬 5.等壓模擬補足蛋白質 X-Ray 結晶失去的氫原子，以及將小分子和水分子加入系統中，水分子以 TIP3PBOX 建立 12Å 此表示蛋白質最外圍的原子距離水分子為 12Å ，經過補氫以及加水後會使得結構產生不穩定的 van der waals 以及 elecctrostatic 作用力，所以我們必須要先進行能量極小化經由 LEaP 處理完的結構，進而才能夠得到較為 local minima 的座標檔案。要進行模擬的初始結構是用 AMBER 程式集中之 LEaP 來產生 prmtop 以及 inpcrd 檔案，其中 prmtop 檔案包含原子電荷、重量；inpcrd 檔案含有原子位置的軌跡檔案。兩者檔案亦包含 3 個 counterion (Na+離子)。進行能量極小化時，首先使用 500 步最陡下降法接著用 500 步共軛梯度法，共計為 1000 步極小化，以期能找到能量相對低點(local minima)。過程中使用 NTB=1 (constant volume periodic boundaries) 。升溫過程當中使用等體積 NTB=1(constant volume periodic boundaries)，與忽略原子 9 Å 以上的作用力(CUT=9)， TI1、TI3 模擬進行在解牛頓方程式時 dt 設定為 2fs，為了避免氫原子震盪頻率過大，使用 SHAKE 對其加入限制(NTC=2，NTF=2)； TI2 則並未使用 SHAKE 對其加入限制(NTC=1，NTF=1)，溫度皆由起始 0K 加熱至 300K 並且使用 Langevin dynamics 且碰撞頻率 2.0p. 來控制溫度 (NTT=3,. gamma_ln=2.0)。等體積模擬是為了看模擬過程當中是否有氣泡的產生。接著使用等壓模擬 constant pressure periodic boundaries 來控制在一大氣壓下的壓力模擬(NTB=2)，忽略原子間 9 Å 外之作用力(CUT=9)，加入 SHAKE，將氫原子與重原子的距離控制在一定的範圍以避免震盪頻率過大。. 18.

(28) 分子動力學模擬過程中，我們需要檢查模擬過程中所得到的結果是否合理，所設的參數條件是否正確，我們會以 2ns 分子動力學模擬結果來分析模擬環境。以下分成四個部分:1.總能量(Etot = potential + kinetic energy)2.溫度 3.密度 4.體積固定體積由 0K 加熱至 300K 的時間內會使得總能量上升，之後溫度壓力固定穩定時，總能量會下降，並且保持穩定的狀態，當從等體積模擬轉變成等壓模擬之時，分子會向外擴散，造成壓力為負值，壓力控制方法為採用控制盒子的體積來形成壓力，但此方法會造成壓力有很大的起伏，不過平均壓力會在 1atm 左右。壓力的控制是用體積的方式，因此在體積的部分也會保持穩定。最後系統密度在水溶液的環境模擬，水分子的密度為 1g/c. ,，若加入蛋白質會使得略微增加。由. 此可見此次分子動力學屬於合理。24. 2.2.2. RESP (Restrained electrostatic potential). 配體上的原子電荷計算上本次使用的方法為 RESP，相較於 AM1-BCC 方法而言，RESP 方法的精準度較為高。因為 AM1-BCC 方法含有半經驗式的參數，假設性的參數使得算出的電荷數值精準度較為低。 RESP 計算電荷是運用古典庫倫模型對應到量子力學分子靜電位能來得到電荷數值，較大之分子欲模擬計算，此方法時常被使用。使用靜電位能的模擬方法電荷通常會重新產生與氫鍵結的錯合物，雖然說事實上這些得到的電荷為：高估了這些交互作用的強度所得之數值，不過將部分的電荷對應到靜電位能時使用 restraint function 能夠使得高估的 ESP-charge 數值衰減。所以此種方法被稱為 RESP-charge。28 使用 RESP 計算電荷會運用到 Gaussian 軟體，此次計算所使用到的版本為 Gaussian 03。Gaussian 為量子化學計算軟體，它以量子化學的原理以數值方法來預測化學分子的性質。當我們使用 Gaussian 軟體時只要提供欲計算的分子三度空間結構，電荷資訊以及理論計算的方法經由大量計算後則可得出電子波函數以及能量並且可以分析其他化學性質。計算所使用的方法以及基底函數(basis set)需相稱可達到平衡，若使用的計算方法與基底函數差異過大，計算所得之數值可能也會導致結果偏離真實的現象。本次計算方法使用 Hartree-Fork 搭配的基底函數為: 6-31G*。Hartree-Fork method 又被稱為 self-consistent field method(SCF)為應用變分 19.

(29) 法來計算多電子體系下波函數方程，變分法是用來處理泛函(Functional)，而得到極值函數。HF 方法通常被用來解決多電子原子或分子下之 time-independent 薛丁格方程式，而去近似基態的波函數以及能量。.. 2.2.3. Functional group transformation. 進行 TI 計算程序欲得到相對自由能，其中一步會進行官能基的變換，也就是說配體會從λ=0 經過變換成為λ=1，必須將原本欲消失的原子變成要得到的原子。產生新的原子接在原來的主體架構上勢必要為合理的構型，符合化學鍵長、鍵角、扭角…及其他化學條件，若直接在 pdb 檔案當中刪除消失原子以及增加生成原子，其實很難建構成為較為合理的結構，此時如能用專業的軟體那麼就能夠較容易的建構出準確分子結構。除此之外尚需考量生成的λ=1 分子，它與λ=0 相同原子部分的位置是否相等?若沒有相等，最後產生之λ=0 以及λ=1 coordinate 檔案當中相同原子部分的各個原子就無法重合，無法符合 TI 計算上的要求。因此經過本實驗室的探索我們採用 DS ViewerPro 6.0 來進行以上的工作，經由此軟體的配合我們可以得到重疊的構型，進而能夠進行接下來的計算。. 2.2.4. Cheng-Prusoff equation. 從下頁表格當中可以看見許多的 ERK 小分子化合物(抑制劑)Ki 值，一般來說藥物化學的論文中有提供小分子化合物 I. 值的資料 29，並且利用 I. 斷小分子化合物對於某酵素抑制效果好壞與否的認定，I. 值來評. 若做出來所得到的數值. 為μM(micromolar)的層次就算不錯，可以達到 nM(nanomolar)就更好，當然藥物化學家不會單就 I. 值來判斷此小分子化合物是否有價值去做更進一步的實驗，. 還要考慮很多其他的因素，不過 I. 是一項指標，值得去參考。. 但是拿在不同實驗室所做出來的 I 不同實驗室拿小分子去測 I. 值資料去做比較是否合理還有待確認！. 值的程序會影響到最後 I. 值的結果，若拿不同實. 驗室的結果直接去做比較那麼可能會有一些誤差的產生。此時如能使用 Ki 值來做比較的話，相較於 I. 會較為理想；因為 Ki 值並不會有 I. 有上述之缺點，所以. 在各個資料庫所蒐集之各個小分子化合物抑制效果比較上 Ki 值相較於 I 如能得到 I. 值皆為從同一實驗室所做出的結果，也是理想的狀態。 20. 來的佳。.

(30) 本篇論文中計算出來得到的數值為相對自由能(ΔΔG) ，需經由公式的轉換才能去計算出論文當中設計的小分子化合物 Ki 數值，以下為公式的推演：. ※{S}代表 substrate、Km 代表平衡常數. ΔG= -RTln. ΔΔ. ΔΔ. Δ. =-RTln. Δ. =-RTln. Δ. =Δ. = -RTln. =Δ. - (-RTln. -Δ. )= Δ. -Δ. ΔΔ. = -RTln. = Δ. -Δ. ΔΔ. = RTln. = Δ. -Δ. ※R=8.314 J. (氣體常數) T=300 K (絕對溫度). 故由以上的推導，可以證明由計算 complex 以及 solution 的相對自由能，並且以下頁表格所提供實驗之起始化合物 Ki 值，能夠得到自行設計化合物的 Ki 值，我們就能進一步跟其他有實驗 Ki 值的化合物進行比較。在此要注意的是:AMBER 系統計算所得到的結果單位為 kcal/mol，單位上要進行換算使得準確數值。. 21.

(31) 圖表 2-410a. 各種類似化合物之活性. 22.

(32) 2.2.5. TTA v.s. STA. 進行熱力學積分得到相對自由能有兩種方式:第一種方式為：TTA (Three Transformation Approach) ；第二種方式為： STA(Single Transformation Approach)。 TTA 的策略為: 運算時，分成三個步驟來進行。 1.將欲進行變換之原子電荷移除。 2.進行凡得瓦變換：移除消失的原子，並且加入新生成的原子。過程當中唯一會使用 softcore potential 之步驟。 3.補上加入原子之電荷。目前在 AMBER 網站上 TUTORIAL 連結當中有一篇介紹 TI 的範例：即採用 TTA 的方法來進行計算，此法相較於 STA 方法來說，它在 AMBER 討論區有較多探討文章，也暗示著兩者相較之下較多人使用 TTA。此方法在第二步使用的 softcore potential 為經過修改之 Lennard-Jones potential30，此過程並不牽涉到電荷的計算，主要是描述凡得瓦作用的能量。 STA 的策略為: 運算時，僅用一步完成計算。也因為使用一步便完成電荷以及凡得瓦作用力之計算，故 softcore potential 與 TTA 之不同。修改過後如同以下: =(1- ) = 此計算的公式是有經過修改的結果，不同於 TTA 方法首先移除原子上的電荷，緊接著進行 van der waals(凡得瓦)變換得到 =1 之分子，最後再將電荷補上新增的原子上。所以計算上的公式也經過略為修改，它所使用的 softcore potential 是 Coulomb potential 來進行能量的計算。不過目前並不確定使用 STA 方法收斂是否相較於分成多個步驟的 TTA 來的快，相關文獻上面並沒有提及到，可能需要再經過更多的研究才能夠得知。不過計算程序上的操作以 STA 來說：它相較於 TTA 是較為的簡便，只需要計算一次即可，綜合以上 STA 還是值得去探討。. 23.

(33) 2.2.6. Comparison of different λ points. 本實驗所計算的 λ 點(積分點)之數目基礎為 λ=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9。一共是九個點。計算 TI 採用 λ 點的數目以現況來說並沒有一個絕對；有的文獻上做的 λ 點比較多，多則達到數十個點；也有做非常少的 λ 點：僅僅只有一個 λ 點。所以 λ 點的選擇五花八門，且不同的研究標的可能情況也不甚相同，也因此目前參照現有的文獻無法下一個定論。採用不同數目的 λ 點，關鍵還是在於我們從 λ=0 至 λ=1 這段過程曲線的平滑程度，若我們在積分的時候各個點所連接成之曲線很平滑，假設線是平的，那麼我們利用很少的點就可以得到較為精確的數值了；相反的若從 λ=0 至 λ=1 的過程各個點所連接成之曲線是非常曲折，那麼 λ 點就需要比較多才能得到較為準確的數值。但是若我們無法預知這段模擬過程當中究竟曲線會是平滑還是曲折，λ 點的數目可採取較為保守的方式從 λ=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9。一共採用九個點。也因此本研究採取的 λ 點數目為九個點。而後再擷取其中幾個點做積分來比較取較少 λ 數目作積分所得之數值是否貼近？還是相距甚遠？. 24.

(34) Chapter 3 3.1. 結果與討論 (Results and Discussions). TTA (Three transformation approach) relative binding free energy( λ points = 9). 3.1.1. TTA：CFO to 82A (solution phase). 圖表 3-1. λ=0 CFO 進行變換形成 λ=1 82A. 首先進行的計算方式為 TTA，將配體 CFO 以及 82A 加上水盒子，計算相對自由能，第一步進行移除電荷：移除紫色圈起來的幾個原子:CF3 以及 OH 官能基。. 圖表 3-2. 配體 λ=0 未含水結構:CFO. 圖表 3-3. 25. 配體 λ=0 含水結構:CFO.

(35) 第二步進行 vdw transformation，將 CF3 轉變成為 F；OH 轉變成為 Cl。. 圖表 3-4. 配體 λ=1 未含水結構:82A. 圖表 3-5. 配體 λ=1 之含水結構:82A. 第三步進行加入電荷的動作，將 F 以及 Cl 的電荷補回。這三部的計算都有產生自由能，而我們再將這三個步驟的自由能數值相加，而後就可得到在水溶液下進行變換的自由能。表格 3-1. solution phase Total Gibbs energy. Kcal/mol. 2eq_prod(0.2ns). 3pro1(0.6ns). TI1. 6.435295. 6.51921. TI2. -6.13642. -4.76397. TI3. -5.329915. -5.3433. Total Gibbs energy. -5.03104. -3.58806. 進行完能量極小化後接續在 10ps 進行等壓升溫到達 300K，不直接使用 300K 而是從 0K 升溫至 300K 的原因是以防突然偏離實驗結構太遠的情況，接續為 200ps 的 equilibration run 以得到 relax 的結構，而後再進行 600ps production run。接著檢察這段模擬過程當中能量是否有得到收斂的情況，與 complex phase 相較之下 solution phase 能量收斂速度較為快速。模擬 λ=0 的狀態下：以 CFO 於水溶液中來說，一共有 4405 個原子；λ=1：82A 於水溶液中則有 4401 個原子，solution phase 下以原子的數量以及模擬系統的複雜度來說，相較於 complex phase 都較為的簡易。所以於這樣的計算中，只需要用到一個亂數就可以了。Solution phase 在整體 TI 的計算量來說，所佔的比例也是相當的低，大概只有百分之五。 26.

(36) 我們看到表格 3-1 的結果可以發現 200ps 時，所計算得到的 total Gibbs energy 為-5.03104kcal/mol ，而 600ps 計算所得之能量為-3.58806kcal/mol，更細膩的去看各個步驟的能量我們可以發現 TI1 以及 TI3 的部分：200ps 以及 600ps 兩者積分出來的數值差異都不大；以 TI1 來說差異約為 0.8kcal/mol、TI3 約為 0.01kcal/mol。可是 TI2 計算的部分差異就較為大了些，兩者差異約為 1.33kcal/mol。原因是第二步為進行 vdw 變換，因此差異會相較 TI1、TI3 來的大。. 3.1.2. TTA：CFO to 82A (complex phase). 圖表 3-6. 2OJJ 緞帶形狀結構紅色部分為 Helix 部分紫色為 strand 部分此圖為經由 UCSF Chimera 所得到之構型圖。. 27.

(37) 圖表 3-7. 2OJJ 加入水盒子後之構型。. 加入水盒子後，系統內共含有五萬多顆原子，而含有 82A 配體的檔案與含有 CFO 配體的檔案相差了 4 個原子。而下圖是觀察配體與周遭胺基酸作用為疏水性或者是親水性的情況，而用黃色圈起來的地方即進行官能基變換的位置。. 圖表 3-8. CFO 以及 82A 在 complex phase 下重疊示意圖. 28.

(38) 而後如同在 solution phase 的計算程序，我們同樣的在 complex phase 下使用同樣的參數條件進行計算，使用 λ 點的數目為 9，且分為 3 個步驟(step)計算。表格 3-2. complex phase TI1&TI3 free energy. Kca/mol. Pro1(0-0.6ns). Pro2(0.6-1.2ns). TI1. 10.45399. 10.33826. TI3. -9.0911. -9.270855. TI1 以及 TI3 計算中我們 mask 的原子分別為：TI1-CF3、OH；TI3-F、Cl。也就是說在 TI1 進行模擬我們是移除掉 CF3、OH 上的電荷，TI3 進行模擬我們是補上 F、Cl 的電荷。我們可以從上表得到 TI1 production run 0-0.6ns(600ps)，所得到的能量為 10.45399kcal/mol ， production run 0.6ns-1.2ns 所得到的能量為 10.33826kcal/mol。進行 TI1 能量部分並未牽涉到 van der waals 轉換，它能量僅用一條軌跡即得到收斂。TI3 的部分也可以發現與 TI1 的情況類似，得到的數值約為 -9 kcal/mol 左右。表格 3-3. complex phase TI 2 使用 6 個亂數平均而得之能量. Random. Pro1(0-0.6ns). Pro2(0.6ns-1.2ns). number(kcal/mol) TI2_71277. -3.496165. -3.44248. TI2_50111. -4.69406. -4.08. TI2_71379. -3.72782. -3.391075. TI2_4113. -4.230975. -3.400715. TI2_275467. -4.193965. -3.49153. TI2_10529. -4.645565. -3.74432. Six Avg. -4.164758. -3.591687. TI2 進行 van dew waals transformation，在這個步驟當中，我們在含有 CFO 之 complex 轉換成含有 82A 之 complex。此步驟做計算所消耗的計算資源很龐大。我們使用了 6 個亂數也就是說 TI2 用不同的亂數計算了 6 次才得到以上的 TI2 平均值。 29.

(39) 亂數的選擇一般使用質數。亂數會產生波茲曼分布，以得到初始速度、初始位置、以及位能函數。用不同的亂數得到不同的結果再求取平均值，相較於使用一個亂數進行很長的一段模擬時間是較為理想。不同亂數之下有時候能量差異較大的情況下會差到 2kcal/mol 左右，這是還可以接受的範圍，當然也有的時候不同的亂數條件下所得到的數值很接近。以表格 3-3 來說：我選用了 6 個亂數，分別為：71277、 50111、71379、4113、275467、10529。Production run 0-0.6ns，差異最大的兩個亂數為 71277、50111，兩者能量差距達到 1.20241 kcal/mol。反之差距最小的兩個亂數為 50111 與 10529，兩者差距為 0.048495 kcal/mol。而為了得到更為收斂的數值，每個亂數都同樣的進行 2 段 0.6ns 的模擬，一共進行 1.2ns。到了第二段 production run 0.6ns-1.2ns 可以發現 6 個亂數下所得到的數值差距有變小的趨勢，最大的差距為 50111 以及 71379 的差值：0.688925 kcal/mol；最小的差距為 71379、4113 之差值：0.00964 kcal/mol。從上可得知，差距明顯的相較於第一段是來的小。由此看來以能量的觀點是有得到良好的收斂。. 30.

(40) 圖表 3-9. 圖表 3-10. 82A 與周圍胺基酸交互作用情形. CFO 與周圍環境胺基酸交互作用情形 31.

(41) 3.1.3. Summary of CFO to 82A 表格 3-4. TI solution & complex phase of total free energy. (kcal/mol). TI1. TI2. TI3. total. solution. 6.51921. -4.76397. -5.3433. -3.58806. complex. 10.33826. -3.591687. -9.270855. -2.524282. 綜合以上的數據 solution phase 我們是採用 0-0.6ns production run 之結果，而 complex phase 則採用 0.6-1.2ns 收斂的結果。Complex phase 下 TI1 以及 TI3 計算時只有使用一個亂數，故直接使用 0.6-1.2ns 的數值；TI2 計算需要經過漫長的收斂，所取的數值為 6 個亂數計算下所得到的平均值。Solution phase 之 TI1、TI2、TI3 三者相加後所得到之數值為-3.58806 kcal/mol；而 complex phase 之 TI1、TI2、TI3 三者相加後所得之數值為-2.524282。按照 Thermodynamic cycle 概念我們可以得到:. ΔΔ. Δ. =Δ. = Δ. -Δ. Delta-delta G= Delta G for complex - Delta G for solution= -2.524282-(-3.58806)=1.063778 kcal/mol 以上的計算可以看到相對自由能為 1.063778kcal/mol，這表示相較於 82A，我們所設計出的小分子化合物 CFO 對於晶體結構 2OJJ 結合能力更為的強。依照計算結果之能量：理論上 CFO 的抑制效果會比 82A 來的好。82A 以實驗對 ERK2 所測出之 Ki 值為 0.002μM=2nM。我們所設計出的 CFO 經由轉換所得到的 Ki 值相較於 82A 會來的更為理想。. 32.

(42) 3.2. Comparison of TTA with STA. 3.2.1. Results of STA (CFO to 82A). 上節 3-1 使用 TTA 計算得到最後的相對自由能。而本節希望能夠簡化計算的程序而使用 STA 來得到相對自由能值。目前探討 TTA 與 STA 的論文不多，文獻上有關 STA 與 TTA 計算之比較不確定是否同樣適合用在本篇的標的。參照 Amber 使用手冊，將計算程序改為 STA，scbeta 參數值使用系統的預設值：12。計算使用的 potential 為 coulomb potential 也就是電荷和凡得瓦計算在同一步驟完成，相異於 TTA 所使用的 modified Lennard-Jones potential。本次得到 STA 相對自由能為計算 complex phase (production run 600ps)以及 solution phase (solution phase 600ps)下之差值。而最終計算出來的數值為：. 表格 3-5. STA(production run 600ps)計算結果. Kcal/mol. STA total. Solution phase. -0.405845. Complex phase. 2.4701. 33.

(43) Delta-delta G= Delta G for complex - Delta G for solution= 2.4701 -(-0.405845)=2.875945 kcal/mol STA 相對自由能值為 2.875945kcal/mol，相較於 TTA(lambda points =9)之相對自由能 1.063778 kcal/mol ，兩者之差值約為 1.812kcal/mol。自由能為 state function(狀態函數)，即無論反應過程路徑是如何變化，所得到的數值為末狀態減掉初狀態所得到數值都會是相同的，以此看來 STA 與 TTA 應該得到的數值要一樣，可是單就 solution phase 來看 STA：-0.405845kcal/mol 與 TTA： -3.58806kcal/mol 相差了 3.182215 kcal/mol，complex phase STA：2.4701kcal/mol 與 TTA：-2.524282kcal/mol 相差 4.994382kca/mol。但若直接看相對自由能 STA 與 TTA 相差 1.812kcal/mol。參考文獻上的積分點數目來增加此次計算的積分點以探討是否能明顯改善。如同上頁所示：將 solution phase lambda 點的數目由原本 9 個點增加為 19 個點，間距由原先 0.1 縮小為 0.05，計算當中增加了：0.05、0.15、0.25、0.35、0.45、0.55、 0.65、0.75、0.85、0.95，一共 10 個點，以此來看是否能夠有效減少與 TTA 相差的數值。最終計算出來得到 solution phase 19 個點的數值約為：0.06597 kcal/mol。與原先 9 個點相差 0.471815kcal/mol。. STA_pro_sol(0-600ps). 25. dV/dλ(kcal/mol). 20 15 10 19lambda. 5. 9lambda. 0 0. 0.2. 0.4. 0.6. 0.8. 1. -5 -10. 圖表 3-11. solution phase 9 個積分點與 19 個積分點 34.

(44) dV/dλ(kcal/mol). STA_pro1_(0-600ps) 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. STA_pro 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9. 圖表 3-12. STA complex phase lambda points=9. Delta-delta G= Delta G for complex - Delta G for solution= 2.4701 -(0.06597)=2.404313 kcal/mol 由以上 19 個點計算出來所得到相對自由能差值與 TTA 差距的確也縮小了(參照下頁表格 3-6)。單就 solution phase 與 TTA solution phase 去做比較也縮小了但是得到的 solution phase 自由能 0.06597 與 TTA solution phase 自由能：. dV/dλ(kcal/mol). -3.58806kcal/mol 還是有段差距。. 50 45 40 35 30 25 20 15 10 5 0 -5 0 -10 -15 -20. 圖表 3-13. comparison of solution. 1_step_sc12_equi 1_step_sc7_equi 3_step 1_step_19_lambda 1_step_sc12_prod 0.2. 0.4. 0.6. 0.8. 1. 1_step_sc7_prod. solution phase 下 TTA 與 STA 及調整參數和不同 λ 數目比較. 依照上圖所示：更改 scbeta 參數由 12 改為參考文獻第 22 篇當中所測試出之優化數值 scbeta 值改為 7。因為按照文獻上所述：我們按照原始 Amber Package m=2， 35.

(45) 所得之數值 2.5 到 3.5 之間，也就是說. 在 6.25 到 12.25 是理想的數值。本. 次實驗所更改之 scbeta 即由 12 改為 7。. 表格 3-6. 更改 scbeta 數值為 7 後之比較. Step. Total gibbs free energy(kcal/mol). 1_step_sc12_equi. -0.554095. 1_step_sc7_equi. 0.16878. 3_step. -3.58806. 1_step_19_lambda. 0.06597. 1_step_sc12_prod. -0.405845. 1_step_sc7_prod. -0.069535. 可發現到無論按照更改過後 scbeta 數值後所得到的數值差距並不大。相較於 scbeta=12 是有些許的改善，不過也可發現，造成這當中的差距之主要原因，應該不在於 scbeta 數值之調整上。而另外在 Amber 手冊當中也有提到：STA 收斂速度相較於 TTA 不一定比較快。故是否有可能因為計算上 STA 所得之結果並未收斂? 值得思考。. 最後經過思考發現在計算程序上似乎在 STA 處理電荷部分是有些不同的，計算 TI 時計算電荷使用到 particle mesh Ewald31，它是在處理 long range electrostatic (長距離靜電)，它有用到 Ewald summation。先前在進行 TTA 時，我們第一步會先將欲變換掉之原子電荷移除，到了第三步再將欲補上之原子的電荷補上；更細微的看待這個過程：第一步移除那些原子電荷後，系統整體會變成非中性，系統為了要將整體變為中性，則會加上失去那些電荷之數值平均補回整體之系統，使之變為中性再進行下一步之 van der waals transformation。而進行 STA 則不同。 36.

(46) 也因此我們嘗試著將每個原子之電荷數值都減半，也就是說減低電荷的效應，看看結果是否如同我們之假設。以下我們進行了電荷減半之 TTA solution phase 以及 STA solution phase，可以發現到與原先並未減半前是有不同的。. 表格 3-7. TTA(solution phase)電荷減半所得之自由能. TTAhalfcharge(kcal/mol). TI1. TI2. TI3. total. Solution_Equi_200ps. 1.692525. -1.477875. -1.91986. -1.70521. Solution_Prod_600ps. 1.720925. -1.41863. -1.928735. -1.62644. 表格 3-8. STA(solution phase)電荷減半所得之自由能. STAhalfcharge (kcal/mol). onestep. Solution_Equi_200ps. -0.006145. Solution_Prod_600ps. 0.009015. 37.

(47) 表格 3-9. TTA(solution phase)電荷減半與未減半之比較. TTA(kcal/mol). Prod_600ps(total). Solution_halfcharge. -1.62644. Solution_normalcharge. -3.58806. solution phase 電荷按照原始來看 STA：-0.405845kcal/mol 與 TTA： -3.58806kcal/mol 相差了 3.182215 kcal/mol。 Solution phase 電荷減半過後 STA：0.009015kcal/mol 與 TTA：-1.62644kcal/mol 相差了約 1.6273415kcal/mol。. 表格 3-10. STA(solution phase)電荷減半與未減半之比較. STA(kcal/mol). Prod_600ps. Solution_halfcharge. 0.009015. Solution_normalcharge. -0.405845. 原先在 solution 電荷尚未減半的例子當中: production run 所相差的數值為:3.182215kcal/mol 而在 solution 電荷已經減半的例子當中:production run 所相差的數值為:1.6273415lcal/mol。發現電荷減半的例子當中 TTA 與 STA 計算數值得到的差距有縮小的跡象與尚未減半的例子當中做比較:相差的數值為:1.5548735kcal/mol。因此證明我們的論點：計算無論是於 solution phase 或者是 complex phase 自由能在 STA 或者是 TTA 下，有差距的主要原因是來自於電荷的因素。. 38.

(48) 3.3. Result of different λ points in TTA. 3.3.1. λ points equal to 9. 計算 TTA 於 complex phase 之下，進行第二步 van der waals 變換時，必須要使用多個亂數來進行模擬才可以得到較為收斂的數值。以下即為進行 TI2 production run (0-0.6ns)的曲線圖：橫軸代表 λ 點，縱軸代表每個 λ 點之自由能。由此圖可看出不同的亂數有些高低起伏的趨勢是相同的，但有些就不太相同了。以下圖來說亂數 275647 相較於其他五條就有點不同，不過我們最後要看的還是積分出來的數. dV/dλ(kcal/mol). 值。. 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. TI2(0-0.6ns). 71277 71379 50111 4113 10529 0.1. 0.2. 圖表 3-14. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. TI2 complex phase(0-0.6ns) λ points=9. 39. 275647.

(49) dV/dλ(kcal/mol). TI2(0.6ns-1.2ns) 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. 71277 71379 50111 4113 10529 0.1. 0.2. 圖表 3-15. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 275647. TI2 complex phases(o.6-1.2ns) λ points=9. 延續 TI2 complex phase (0-0.6ns) 的結果，以此可以發現 TI2 complex phase(0.6-1.2ns)的結果曲線較為重疊，且最後積分出來的結果 6 個亂數值上的差異. dV/dλ(kcal/mol). 也較為的小，即表示能量值趨向於收斂。. 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. complex phase TI1,TI3. TI1(0.6-1.2ns) TI3(0.6-1.2ns). 0.1. 0.2. 0.3. 圖表 3-16. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. complex phase(0.6-1.2ns) λ points=9. TI1 complex phase(0.6-1.2ns)，圖形當中可發現曲線相較於 TI2，以及 TI3 都更為的陡峭，以能量的觀點來說，TI1 以 1 個亂數去模擬至 production run(0.6-1.2ns) 是已經收斂了，但是之後若取較少的 λ 點數目來做積分，在 λ 點在 0.1 至 0.3 的地 40.

(50) 方可能會造成最後積分結果相較於 λ 點數為 9 的數值差異較大。也因此 λ 為 0.1 到 0.3 的區間，要更為留意。. TI3 complex phase(0.6-1.2ns)中可觀察到相較於 TI1 complex phase(0.6-1.2ns) 曲線是平滑的許多，是 TI 裡面 3 個步驟當中最為平滑的曲線。故若取較少的 λ 點. dV/dλ(kcal/mol). 對於 TI3 積分情況來講，所造成積分出來的數值差異會較為小。. 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. solution phase. 3pro1_TI1 3pro1_TI2 3pro1_TI3 0.1. 0.2. 圖表 3-17. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. solution phase (0-0.6ns) λ points=9. TI2 solution phase (0-0.6ns)當中可發現曲線是很平滑的，也僅僅在 λ 點為 0.7 開始到 0.9 呈現平滑上升，其餘從 λ=0.1 至 λ=0.6 點幾乎是在-5kcal/mol 上下游走。. TI1 solution phase (0-0.6ns)當中可發現曲線相較於 TI1 complex phase (0.6-1.2ns) 來說趨勢是一致向下的。並無呈現高低起伏的狀態，故以此推論在減少 λ 點時，對於最後積分出來的數值影響應該不大。. TI3 solution phase (0-0.6ns)也可發現如同 TI1 般，整體的趨勢也是向下，也沒有高低起伏的情況，故可推論減少 λ 點對於積分出來的數值影響也不大。. 41.

(51) 3.3.2. λ points equal to 5 表格 3-11. λ points=5 能量總整理. Kcal/mol. TI1. TI2. TI3. total. solution. 6.5119375. -4.721535. -5.3214125. -3.53101. complex. 10.19492. -3.5085. -9.28608. -2.59966. Delta-delta G= Delta G for complex - Delta G for solution= -2.59966-(-3.53101)=0.93135 kcal/mol 最後從 λ 點數為 5 的結果當中我們可以看到在 complex TI1 與 9 個點的結果相差約 0.14kcal/mol，是最多的。λ 點數為 5 所取的積分點 0.1,0.3,0.5,0.7,0.9，依據圖表 3-16 的圖形來說從 λ 點 0.1 至 0.3 之間，0.2 為突起的點，因為所取的點沒有 0.2，所以也導致少積分了一小塊面積，所以相差的數值與其他各組比較是最. dV/dλ(kcal/mol). 多的。. 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. TI2(0-0.6ns) 71277 71379 50111 4113 10529 0.1. 圖表 3-18. 0.3. 0.5. 0.7. 0.9. 275647. TI2 complex phase(0-0.6ns) λ points=5. TI2 complex phase(0-0.6ns)沒有取 λ=0.2 這個點，可發現到 0.1 到 0.3 是平的，但是因為 TI2 λ =0.1 到 0.3 沒有明顯的高低起伏故還不夠成巨大影響，其它的地方都不太偏離原先 λ 點為 9 的趨勢。那麼積分出來與 9 個點也就會比較接近。. 42.

(52) dV/dλ(kcal/mol). 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. TI2(0.6ns-1.2ns) 71277 71379 50111 4113 0.1. 0.3. 圖表 3-19. 0.5. 0.7. 10529. 0.9. 275647. TI2 complex phase(0.6-1.2ns) λ points=5. TI2 complex phase(0.6-1.2ns)與前面圖表 3-18 做比較差異性不大同樣也只有亂數為 275647 與其他條的曲線較為不同，但是積分出來的結果與其他亂數值差異性. dV/dλ(kcal/mol). 不大，並非 6 個亂數中得到積分之最大值或者是最小值。. 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. complex phase TI1,TI3. TI1(0.6-1.2ns) TI3(0.6-1.2ns) 0.1. 圖表 3-20. 0.3. 0.5. 0.7. 0.9. complex phase(0.6-1.2ns) λ points=5. TI1 complex phase(0.6-1.2ns)少了 λ=0.2 這個點後，原先要先往上再往下，變成直接向下，表示這並沒有準確的描述出 λ 點 0.1-0.3 之間模擬過程當中曲線上下的變化，也因為如此使得積分出來的數值會差異較大。 43.

(53) TI3 complex phase(0.6-1.2ns) ，λ 點數改成 5 以後與 λ 點數為 9 所呈現的趨勢. dV/dλ(kcal/mol). 是一致的，也就是說利用這五個點就可以得到滿接近 9 個點所積分出來的數值了。. 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. solution phase. 3pro1_TI1 3pro1_TI2 3pro1_TI3 0.1. 0.3. 圖表 3-21. 0.5. 0.7. 0.9. solution phase (0-0.6ns) λ points=5. solution phase 下所觀察到的 λ 各點變化都很小，且趨勢都一致。以 solution phase TI2 來說，它跟 complex phase TI2 下的情況差異頗大，solution phase TI2 曲線變化較少，故於 solution TI2 減少 λ 點的數目對於最後的積分數值也同樣的不至於有太大影響。. TI1 solution phase λ 點數為 5 也同樣因為與 λ 點數為 9 的 solution phase 情況類似，所以並不會造成積分數值有較為偏差的情況。. TI3 solution phase λ 點為 5 也與 λ 點為 9 的 solution phase 趨勢幾乎一致，所以最後積分出來的數值與 λ 點為 9 的數值頗為接近。. 44.

(54) 3.3.3. λ points equal to 4 表格 3-12. λ points=4 能量總整理. Kcal/mol. TI1. TI2. TI3. total. solution. 6.31518. -4.69984. -5.2722. -3.65686. complex. 11.72433. -3.6823. -9.01265. -0.97057. Delta-delta G= Delta G for complex - Delta G for solution= -0.97057-(-3.65686)=2.68629 kcal/mol 於 λ 點數為 4 的情況下，再度發現 TI1 complex phase 與 λ 點數為 9 時差異明顯。λ 點數為 4 來看 TI1complex phase 所積分出來的數值約為 11.72 kcal/mol，與 λ 點數為 9 的積分數值 10.33 差異較大，相差大約 1.4kcal/mol。本次 λ 點數為 4 所取的各點為：0.2、0.4、0.6、0.8，TI1complex phase 中，λ=0.1 與 λ=0.3 相較於 λ=0.2 積分點的數值是比較小的，若只取 λ=0.2 那勢必會造成積分出來的數值相較於 λ 點數為 9 大上許多，所以最終造成與 λ 點數為 9 差異達到約 1.4kcal/mol。平均而言也由於 λ 點數的減少，所以使得相較於 λ 點數為 5 來說，complex phase & solution phase 之 TI1、TI2、TI3 與 λ 點數為 9 差異都來的大，也因此最終加總出來的相對自由能值就比 λ 點數為 9 差距達到約 1.6kcal/mol 左右。. 45.

(55) dV/dλ(kcal/mol). 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. complex phase TI1,TI3. TI1(0.6-1.2ns) TI3(0.6-1.2ns) 0.2. 圖表 3-22. 0.4. 0.6. 0.8. complex phase(0.6-1.2ns) λ points=4. TI1 complex phase 從 λ=0.2 至 λ=0.4 是呈現平滑狀，與 λ 點數為 9 不同，這是造成這邊數值會有比較大差異的主要原因。先前 λ 點數為 9 的情況下 λ=0.1-0.3 之間有呈現倒 V 自行狀態，故以上的線條與先前 λ 點數為 9 不吻合。. TI3 complex phase 與 λ 點數為 9 的曲線所呈現的情況類似，所呈現的曲線平. dV/dλ(kcal/mol). 緩，也因此與 λ 點為 9 的積分數值差異不大。. 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. TI2(0-0.6ns) 71277 71379 50111 4113 0.2. 圖表 3-23. 0.4. 0.6. 0.8. TI2 complex phase(0-0.6ns) λ points=4 46. 10529 275647.

(56) 一般說來 TI2 complex phase 之各個積分點會較為上下起伏，所以減少 λ 點的數目就可能會較為容易有偏差的情形。以上圖來說因為 λ 點數只有 4 個。也就比較會略有誤差多一些的情況，當然還是要看高低起伏上下的情況是否劇烈，若只. dV/dλ(kcal/mol). 是微幅的上下起伏那麼可能誤差會小一點，反之則會形成較大的影響。. 50 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20. TI2(0.6ns-1.2ns) 71277 71379 50111 4113 10529 0.2. 圖表 3-24. 0.4. 0.6. 0.8. 275647. TI2complex phase (0.6ns-1.2ns) λ points=4. TI2complex phase 進行到(0.6-1.2ns)之時，可發現曲線較(0-0.6ns)的平滑，以能量的觀點而言 TI2(0.6-1.2ns)已經較為收斂了，λ 點數為 9：TI2 complex phase 圖型曲線已經較為平滑了，也因此減少點數對於 λ 點數為 4 的影響就會降低。. 47.