以熱力學積分分子動力學模擬計算小分子抑制劑與 GSK-3β 激酶的相對結合自由能

全文

(1)國立臺灣師範大學化學系碩士論文 Department (Graduate institute) of Chemistry, National Taiwan Normal University Master Thesis 指導教授：孫英傑博士 Advisor: Ying-Chieh Sun, Ph.D. 以熱力學積分分子動力學模擬計算小分子抑制劑與 GSK-3β 激酶的相對結合自由能 Computation of Relative Binding Free Energy of GSK-3β Kinase-inhibitor Complexes Using Thermodynamic Integration Molecular Dynamics Simulation 研究生：李興洲 Graduate: Hsing-Chou Lee 中華民國一百零二年六月 June, 2013.

(2) 誌謝師大六年的時光，是人生中轉變、成長的重要階段與回憶。感謝這段時間中給予指導、教誨的師長，以及陪伴我一起走過的同學、朋友和家人。誠摯的感謝指導教授孫英傑博士在研究、學習以及論文撰寫上的指導。帶領我進入原先不熟悉的化學領域，老師悉心的教導使我得以一窺分子模擬領域的深奧，不時的討論、指點研究方向，讓我逐漸的成長、獲益良多。感謝老師在研究態度、人生經驗以及待人處事的指導、分享，是我一生永遠受用的學習經驗。感謝口試委員國立臺灣師範大學生命科學系李冠群教授以及輔仁大學化學系楊小青教授對於論文悉心的指導，使得研究成果能更正確完整的呈現。感謝國家科學委員會的經費以及國家高速網路計算中心在計算資源以及技術上的提供與協助。感謝實驗室的研究夥伴育劭、惠真、嘉仁、玗修、文綺，在研究過程中彼此的討論、激勵，不時的激盪出更深入的想法。感謝實驗室學長姐冠緯、博晉、永富、孟璇、彥樓在碩一階段對於研究以及課業上的指導與經驗分享。最後感謝父母親對於我的栽培與支持，讓我在自由的環境中追尋自己對於學業、人生中的目標。感謝哥哥在學習經驗上的分享。因為您們使得這篇論文能夠順利的完成。. 李興洲. 謹致. 2013 年 6 月. i.

(3) 中文摘要人體腦部神經元中 GSK-3β 激酶蛋白過度磷酸化下游蛋白質 tau 蛋白被認為是引發阿茲海默症可能的病理機制之一。本研究以熱力學積分-分子動力學模擬程序首先計算有實驗結構，及結合能實驗數據的 GSK-3β-小分子抑制劑複合體之相對結合自由能，以檢視計算程序是否再現實驗結果，結果得到誤差為 1.3 kcal/mol。為尋找更好的抑制劑分子，我們從已知的抑制劑分子 ZRM 為基礎，進行官能基轉換，發展 ZRM 之結構相似物並進行的相對結合自由能的計算。最好的計算結果中， ZRM 分子與結構相似物間對於 GSK-3β 之相對結合自由能為 -0.8 kcal/mol，顯示相似物之抑制能力優於 ZRM 分子。我們討論置換官能基與預測的結合自由能的關係。這些結果應有助於在計算模擬以及實驗領域中更進一步發展有效的 GSK-3β 抑制劑。. 關鍵字：分子動力學模擬、熱力學積分、相對結合自由能、肝醣合成酶激酶-3β. ii.

(4) ABSTRACT GSK-3β kinase phosphorylates the downstream protein tau protein. It is thought that, over-phosphorylation of this reaction is one of the pathological mechanisms inducing Alzheimer's disease. In the present study, we calculated the relative binding free energy of GSK-3β-inhibitor complexes using Thermodynamic Integration Molecular Dynamics Simulation and compared it to the experimental result for those complexes of available experimental structures and binding affinity data first in order to investigate how well the computational protocol reproduces experimental data. The computed result was found to deviate from experimental result by 1.3 kcal/mol, showing good agreement. With the computational protocol, we further predict a number of analogs of ZRM ligand in searching for better inhibitors. The best relative binding free energy of GSK-3β-inhibitor complex between ZRM and the analogue is -0.8 kcal/mol, showing better inhibition ability. Substitutions of analogs and their relations with computed relative binding free energy were discussed. These results should be useful for the further development of GSK-3β inhibitors in the areas of computer simulation and experimental assays in the future.. Keywords: Molecular Dynamics Simulation, Thermodynamic Integration, Relative Binding Free Energy, GSK-3β. iii.

(5) 目錄 (CONTENTS) 口試委員會審定書 ........................................................................................................... # 誌謝 ....................................................................................................................................i 中文摘要 .......................................................................................................................... ii ABSTRACT .................................................................................................................... iii 目錄 (CONTENTS) .........................................................................................................iv 圖目錄 (LIST OF FIGURES) ....................................................................................... vii 表目錄 (LIST OF TABLES) ...........................................................................................ix Chapter 1. 緒論............................................................................................................ 1. 1.1. 前言 ................................................................................................................ 1. 1.2. GSK-3 (Glycogen synthase kinase 3) ............................................................. 1. 1.3. 阿茲海默症與 GSK-3β ................................................................................. 2. 1.4. 標靶藥物篩選 ................................................................................................ 2. 1.5. GSK-3β 抑制劑研發現況 .............................................................................. 3. 1.6. 研究目標 ........................................................................................................ 5. Chapter 2. 理論與方法 ............................................................................................... 6. 2.1. 分子動力學模擬 (Molecular Dynamic, MD) ............................................... 6. 2.2. 熱力學積分 (Thermodynamic Integration, TI) ............................................. 6. 2.3. Three-transformation approach, TTA ............................................................. 9. 2.4. 分子模擬系統的建構 .................................................................................. 10. 2.5. 熱力學積分-分子動力學 (thermodynamic integration-molecular dynamic, TI-MD) 的模擬計算 ................................................................................... 13 iv.

(6) 2.5.1. 計算步驟 ............................................................................................. 13. 2.5.2. 計算相關設定 ..................................................................................... 14 結果與討論 ............................................................................................. 16. Chapter 3 3.1. 已知抑制劑分子 ZRL, ZRK 間相對結合自由能計算 ............................... 16. 3.2. ZRL → ZRK 凡得瓦參數調整 ................................................................. 20. 3.3. ZRM, 結構類似物(Mutants)間相對結合自由能計算 ............................... 23. 3.4. 3.3.1. Mu1 ..................................................................................................... 24. 3.3.2. Mu2 ..................................................................................................... 26. 3.3.3. Mu3 ..................................................................................................... 26. 3.3.4. Mu4 ..................................................................................................... 27. 3.3.5. Mu5 ..................................................................................................... 27. 3.3.6. Mu7 ..................................................................................................... 28. 3.3.7. Mu9 ..................................................................................................... 30. 3.3.8. ZRM → M14~M24 ........................................................................... 36. ZRM, M21 與 GSK-3β 複合結構於 1200ps 模擬過程中的軌跡分析。43 3.4.1. 羥基苯基與主結構之間二面角 (dihedral angle) 分析 ................... 43. 3.4.2. M21 分子與周圍原子的氫鍵分析 ................................................... 45. 3.5. M19 與 GSK-3β 複合結構構型分析 ....................................................... 48. 3.6. ZRM, 結構類似物(Mutants)間相對結合自由能計算結果整理 ............... 51 3.6.1. 抑制能力明顯低於 ZRM，ΔΔGbind > 0.5 kcal/mol .......................... 51. 3.6.2. 抑制能力略低於 ZRM，0 < ΔΔGbind ≦0.5 kcal/mol ...................... 51. 3.6.3. 抑制能力略高於 ZRM，-0.5 ≦ ΔΔGbind < 0 kcal/mol ................... 51. 3.6.4. 抑制能力明顯高於 ZRM，ΔΔGbind ≦ -0.5 kcal/mol ..................... 51 v.

(7) 3.7 Chapter 4. 22A, 結構類似物(Mutants)間相對結合自由能計算................................. 55 結論.......................................................................................................... 56. REFERENCES ................................................................................................................ 58 附錄 ................................................................................................................................. 63. vi.

(8) 圖目錄 (LIST OF FIGURES). Fig. 1.1. 已知抑制劑分子 ZRM, ZRL, ZRK .............................................................. 4. Fig. 2.1. 2 相似系統與過渡狀態 ................................................................................. 7. Fig. 2.2. 熱力學積分中 ΔG 的計算 ........................................................................... 7. Fig. 2.3. 熱力學循環(S. Gentile et al, 2011)7 .............................................................. 8. Fig. 2.4. Three-transformation approach, TTA ............................................................. 9. Fig. 2.5. 3ZRM 原始 pdb 構型 ................................................................................ 10. Fig. 2.6. ZRM- GSK-3β (B chain)複合體 .................................................................. 11. Fig. 2.7. ZRM 分子與結晶水、GSK3-β 之間形成氫鍵(G. Gentile et al 2011)8 .... 12. Fig. 2.8. TI-MD 計算流程 ......................................................................................... 14. Fig. 3.1. ZRL → ZRK 熱力學循環 .......................................................................... 16. Fig. 3.2. ZRL → ZRK ΔGcomplex TI2 600ps 計算結果與模擬時間關係圖 ............ 18. Fig. 3.3. ZRL → ZRK ΔGcomplex TI2 600ps 調整凡得瓦參數，計算結果與模擬時間關係圖 ...................................................................................................... 21. Fig. 3.4. ZRM 分子結構與原子編號......................................................................... 23. Fig. 3.5. Mu1 分子結構與周遭 GSK-3β 蛋白質環境 ........................................... 24. Fig. 3.6. ZRM → Mu1 ΔGcomplex TI2 600ps 計算結果與模擬時間關係圖 ........... 25. Fig. 3.7. Mu1 分子結構與周遭 GSK-3β 蛋白質環境 ........................................... 26. Fig. 3.8. Mu3 分子結構與周遭 GSK-3β 蛋白質環境 ........................................... 26. Fig. 3.9. Mu3 分子結構與周遭 GSK-3β 蛋白質環境 ........................................... 27. Fig. 3.10. Mu5 分子結構與周遭 GSK-3β 蛋白質環境 ........................................... 27 vii.

(9) Fig. 3.11. Mu7 分子結構與周遭 GSK-3β 蛋白質環境 ........................................... 28. Fig. 3.12. ZRM → Mu7 ΔGcomplex TI2 600ps 計算結果與模擬時間關係圖 ........... 29. Fig. 3.13. Mu9 分子結構與周遭 GSK-3β 蛋白質環境 ........................................... 30. Fig. 3.14. ZRM → Mu9 ΔGcomplex TI2 600ps 計算結果與模擬時間關係圖 ........... 31. Fig. 3.15. ZRM 分子與 GSK-3β 之結合模式 ............................................................. 33. Fig. 3.16. Mu3, Mu9 與 M24 分子結構.................................................................... 36. Fig. 3.17. M19 與 M21 分子結構 ............................................................................. 37. Fig. 3.18. ZRM → M17 ΔGcomplex TI2 1200ps 計算結果與模擬時間關係圖 ......... 39. Fig. 3.19. ZRM → M19 ΔGcomplex TI2 1200ps 計算結果與模擬時間關係圖 ......... 40. Fig. 3.20. ZRM → M21 ΔGcomplex TI2 1200ps 計算結果與模擬時間關係圖 ......... 41. Fig. 3.21. ZRM → M23 ΔGcomplex TI2 1200ps 計算結果與模擬時間關係圖 ......... 42. Fig. 3.22. M21 分子結構 .............................................................................................. 43. Fig. 3.23. M21 分子結構與二面角示意圖 ................................................................. 44. Fig. 3.24. ZRM, M21 與 GSK-3β 複合結構中，二面角與時間關係圖。 ............. 45. Fig. 3.25. M21 於 1200ps 模擬中與周圍原子形成氫鍵的比例 ................................ 46. Fig. 3.26. 蛋白質結合位深處水分子與周遭原子氫鍵形成的比例 .......................... 47. Fig. 3.27. WAT 6194 與 M21 以及 GLN 150 之間氫鍵形成的比例 ......................... 48. Fig. 3.28. M19 TI3 λ=0.9 1200ps 計算結果氫鍵分析 ............................................... 49. Fig. 3.29. M19 TI3 λ=0.9 200ps 計算結果氫鍵分析 ................................................. 50. viii.

(10) 表目錄 (LIST OF TABLES). Table 3.1. ZRL → ZRK 單一軌跡計算結果 (kcal/mol) .......................................... 17. Table 3.2. ΔGcomplex TI2 不同軌跡計算結果 ............................................................... 17. Table 3.3. ZRL → ZRK ΔΔGbind 600ps 計算結果 ..................................................... 19. Table 3.4. ZRL → ZRK ΔΔGbind 600ps 調整凡得瓦參數，計算結果 ..................... 22. Table 3.5. ZRM → Mu1 ΔΔGbind 600ps 計算結果 .................................................... 25. Table 3.6. ZRM → Mu7 ΔΔGbind 600ps 計算結果 .................................................... 29. Table 3.7. ZRM → Mu9 ΔΔGbind 600ps 計算結果 .................................................... 32. Table 3.8. ZRM → Mutants 計算結果 1 ................................................................... 34. Table 3.9. ZRM → Mutants 計算結果 2 ................................................................... 35. Table 3.10 ZRM → Mutants 計算結果 3 ................................................................... 38 Table 3.11 ZRM → M17 ΔΔGbind 1200ps 計算結果 .................................................. 39 Table 3.12 ZRM → M19 ΔΔGbind 1200ps 計算結果 .................................................. 40 Table 3.13 ZRM → M21 ΔΔGbind 1200ps 計算結果 .................................................. 41 Table 3.14 ZRM → M23 ΔΔGbind 1200ps 計算結果 .................................................. 42 Table 3.15 ΔΔGbind > 0.5 kcal/mol 之 Mutant 分子 ....................................................... 52 Table 3.16 0 < ΔΔGbind ≦0.5 kcal/mol 之 Mutant 分子 ............................................... 53 Table 3.17 -0.5 ≦ ΔΔGbind < 0 kcal/mol 之 Mutant 分子 ........................................... 54 Table 3.18 ΔΔGbind < -0.5 kcal/mol 之 Mutant 分子 ..................................................... 54 Table 3.19 22A → Mutants 計算結果 ........................................................................ 55. ix.

(11) Chapter 1. 緒論. 前言. 1.1. 近年來，標靶藥物的發展為神經退化性疾病 (Neurodegenerative disease)、癌症等疾病的治療帶來了新的曙光。透過對於細胞生物學的了解，針對致病機制中關鍵的訊息傳遞因子或受體蛋白進行藥物設計，從中調節細胞的生理機能，進而達到治療、改善疾病的效果。人類腦部神經元中的 GSK-3β 活性與阿茲海默症患者腦部神經退化具有相當程度的關聯性. 12. 。本研究利用分子模擬 (molecular simulation) 的方法尋找有效抑. 制 GSK-3β 活性的小分子藥物，期望為阿茲海默症治療的研究有所貢獻。. 1.2. GSK-3 (Glycogen synthase kinase 3) GSK-3 Glycogen synthase kinase 3 是一種絲氨酸 / 蘇胺酸蛋白質激酶. (serine/threonine protein kinase) 可催化蛋白質上特定胺基酸的磷酸化反應。學術界最初對於 GSK-3 系列蛋白質功能的了解在於，GSK-3 能磷酸化肝醣合成酶 (Glycogen synthase)，磷酸化會導致肝醣合成酶失去活性 5，藉此機制調控細胞及人體中血糖與肝醣的比例。而後，更廣泛的研究發現 GSK-3 可磷酸化細胞中超過 40 種蛋白質 13, 27，可見 GSK-3 在細胞生理機制中的複雜及多樣性。哺乳類的 GSK-3 編碼於 GSK3α 以及 GSK3β 兩個基因中。近期的研究指出 GSK-3 與多種人體疾病相關，包含第二型糖尿病 (Type II diabetes)、阿茲海默症 (Alzheimer’s Disease)、發炎反應 (Inflammation)、癌症 (Cancer)、躁鬱症 (Bipolar disorder) ，目前我們合作的實驗室特別針對 GSK-3β 的抑制劑運用在阿茲海默症新藥的研發上進行相關研究。 1.

(12) 阿茲海默症與 GSK-3β. 1.3. 阿茲海默症是一種慢性神經退化性疾病，目前全世界約有 2 千 5 百萬人受此疾病困擾。此症主要造成漸進式的認知能力退化，如：記憶、方向感等能力的喪失。阿茲海默症在病理上有二個主要能觀察到的特徵，其一為細胞外類澱粉蛋白的沉積 (amyloid deposits) ，其二為神經纖維糾結 (neurofibrillary tangles, NFT) ，而其中神經纖維糾結主要導因於 tau 蛋白的過度磷酸化 (hyperphosphorylation) 。在人腦的神經元細胞中含有一定數量的 GSK-3β 能對 tau 蛋白進行磷酸化的作用，tau 蛋白與細胞中微管 (microtubule) 的生成相關 29，tau 蛋白的過度磷酸化會致使其與微管分離，分離的 tau 蛋白於神經元中累積與聚集，生成成對的螺旋微絲近一步形成神經纖維糾結，最後導致神經元細胞的死亡 12。 GSK-3β 在神經元細胞中的活性，與阿茲海默症的病理現象具有一定程度的關聯性，使其成為令人感興趣的藥物標的，抑制 GSK-3β 能減緩後續與此症相關的病理現象。GSK-3β 對於 tau 蛋白的磷酸化的過程中，首先 ATP 進入 GSK-3β 的結合區 (binding pocket) 中，藉由 GSK-3β 催化將 ATP 上的磷酸根轉移至 tau 蛋白，生成磷酸化的 tau 蛋白，以及 ADP。大部分的 GSK-3β 抑制劑的作用原理是與 ATP 競爭結合區，使 ATP 無法與 GSK-3β 結合，近而阻斷後續的磷酸化作用。. 1.4. 標靶藥物篩選良好、具治療效用的 GSK-3β 抑制劑分子需考慮許多因素，包含 ADME-Tox28. 即藥物於人體中的吸收 (absorption)、分布 (distribution)、代謝 (metabolism)、排泄 (excretion)、毒性 (toxicity)，而作用於中樞神經系統的藥物，因血腦障蔽 (blood– 2.

(13) brain barrier) 的阻隔需格外考量藥物穿透細胞膜的能力。此外，藥物對於多種蛋白質的選擇性和專一性 (selectivity and specificity) 亦是重要考量因素. 14. ，激酶. (kinase) 這類型的蛋白質主要負責磷酸化下游蛋白，在細胞的訊息傳遞路經中，常藉由磷酸化級聯 (phosphorylation cascade)24 的方式，從細胞膜表面的受體開始，層層的磷酸化下游蛋白質，達成訊息傳遞的目的，藥物分子對於蛋白質的選擇性和專一性不足可能導致多種標靶以外的蛋白質同樣遭受抑制，造成過強的細胞毒性以及其他副作用。根據 Pharmaceutical Research and Manufacturers of America, PhRMA 的資料 22，如(附圖 1)中所示，一顆新藥的研發至最後上市，約需要 10~15 年的時間。而平均研發一顆新藥所需投入的經費則逐年增長，從 1970 年代的 1 億 4 千萬美元，至本世紀初期上升至可觀的 12 億美元(包含研發失敗所投入的成本)，30 年間有將近十倍的增長。PhRMA 成員每年投資在研究發展的花費，在 2011 年估計將高達近 500 億美元 22。由於臨床研究曠日廢時，如何在藥物開發 (Drug discovery) 以及臨床前實驗 (Preclinical trail) 的階段更有效率地從眾多的候選藥物中，發現具有潛力的分子，可能成為縮短藥物研發進程以及所需經費的關鍵。隨著計算機效能的提升以及計算理論的發展，逐漸提升了電腦輔助藥物設計在新藥的研發過程中的重要性。預測未知藥物分子對於特定蛋白質抑制能力的相關計算理論、誤差以及相關應用在近期文獻 9, 21 中被加以探討。整體來說，為減少計算結果的誤差，通常需使用理論精確度 (Theoretical accuracy) 較高的方式，但相對需要較長的計算時間以及較多的計算資源(附圖 2)。. 1.5. GSK-3β 抑制劑研發現況近十年內，GSK-3β 相關的抑制劑在專利以及學術文獻中被大量的加以報導。 3.

(14) 與 GSK-3β 相關的專利申請案件，從 2000 年的 6 件逐年上升至 2008 年的 136 件， GSK-3β 抑制劑的開發及其治療潛力近年逐漸在學界以及製藥產業中受到重視。相關專利中以與阿茲海默症治療相關為最大宗，其次分別為癌症與糖尿病 23。文獻中 23，以 KinaseChemBioBase™搜尋學術文獻以及專利中報導的 GSK-3β 抑制劑分子，高達 45,520 個分子被認為具有抑制能力，但僅 8,193 個分子具有量化的抑制能力資訊，被報導的分子中半數最大抑制濃度(half maximal inhibitory concentration, IC50)值大多落在 10 nM - 1μM 的範圍之間，僅少數約 12 個分子具有較高的抑制能力 IC50 < 1 nM。作者分析抑制劑分子的構型，其中嘧啶(pyrimidine) 相關衍生物占最大的比例，其次則為吡啶(pyridine)相關衍生物。近期文獻中 8，報導了 2-(4 吡啶基)噻吩並吡啶酮(2-(4-pyridyl)thienopyridinone) 相關衍生物對於 GSK-3β 的抑制能力，抑制劑分子的 IC50 值介於 5 nM – 500nM 之間。作者使用 X 射線晶體學(X-ray crystallography)的方式分析結晶構型，探討其中三個抑制劑分子與 GSK-3β 的結合模式(binding mode)。文獻中抑制劑分子與 GSK-3β 複合體的結晶構型於 RCSB Protein Data Bank10 資料庫中的 PDB code 為 3ZRM, 3ZRL, 3ZRK，相對應的抑制劑分子 Ligand id 分別為 ZRM, ZRL, ZRK，分子構型如 Fig. 1.1。. Fig. 1.1 已知抑制劑分子 ZRM, ZRL, ZRK 4.

(15) 1.6. 研究目標本研究期望以文獻中報導的已知 GSK-3β 抑制劑分子 ZRM 為基礎，進行小. 部分官能基的改變，以熱力學積分分子動力學模擬(Thermodynamic Integration Molecular Dynamics Simulation, TI-MD)，預測變動官能基是否有助於提升小分子藥物的抑制能力，為後續的生化實驗提供具有潛力的候選藥物。. 5.

(16) 理論與方法. Chapter 2 2.1. 分子動力學模擬 (Molecular Dynamic, MD) 對於任一系統給予一質點初始座標、速度，並定義周圍的作用力場 (force field). 求解牛頓方程式： 𝐹 = 𝑚𝑟̈ = −. 𝑑𝑉 𝑑𝑟. 可以預測該質點在一段時間中運行的軌跡，解決我們所熟悉的物理習題，如：彈道軌跡、行星運動等。分子動力學模擬則將此一概念應用於微觀的分子系統中，將各個原子視為質點，初始速度則根據模擬溫度以及波茲曼分佈以亂數決定，本研究所使用的分子動力學模擬軟體 Amber113 中對於作用力場的定義為： 𝑉(𝑟) = ∑ 𝐾𝑏 (𝑏 − 𝑏0 )2 + ∑ 𝐾𝜃 (𝜃 − 𝜃0 )2 + 𝑏𝑜𝑛𝑑𝑠. 𝑎𝑛𝑔𝑙𝑒𝑠. ∑. (𝑉𝑛 ⁄2)(1 + 𝑐𝑜𝑠[𝑛𝜙 − 𝛿]). 𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙𝑠. + ∑ (𝐴𝑖𝑗 ⁄𝑟𝑖𝑗12 ) − (𝐵𝑖𝑗 ⁄𝑟𝑖𝑗6 ) + (1⁄4𝜋𝜖0 ∙ 𝑞𝑖 𝑞𝑗 ⁄𝑟𝑖𝑗 ) 𝑛𝑜𝑛𝑏𝑖𝑗. 分別代表鍵的伸縮 (bonds)、鍵角 (angels)、扭角 (dihedrals)、非鍵結作用力 (nonbij 包含凡得瓦作用力以及靜電作用力)。分子模擬系統相當的龐大、複雜，模擬軟體透過演算法，求得數值解，預測系統中各個原子在一段時間中移動的軌跡。. 2.2. 熱力學積分 (Thermodynamic Integration, TI) 本研究使用熱力學積分 15 原理計算二相似系統之間自由能的差值 (ΔG) 𝑉(𝜆) = (1 − 𝜆)𝑉0 + 𝜆𝑉1. λ=0 以及 λ=1 分別代表二相似系統，λ=0 的系統位能為 V0 ，λ=1 則為 V1，λ 介. 6.

(17) 於 0~1 之間則代表 2 個相似系統之間的過渡狀態，如 Fig. 2.1 所示。過渡狀態的位能函數 V(λ) 可由上式求得。. Fig. 2.1 2 相似系統與過渡狀態 1. ∆𝐺 = ∫ ⟨ 0. 𝛿𝑉 ⟩ 𝑑𝜆 𝛿𝜆 𝜆 𝑑𝑉. 二相似系統之間的自由能差異可由上式求得，如果以 λ 為 X 軸，⟨ 𝑑λ ⟩λ 為 Y 軸作圖，則曲線下的面積即為 ΔG. Fig. 2.2 熱力學積分中 ΔG 的計算熱力學積分可應用於藥物分子與蛋白質結合能力 7,. 7. 31. 、蛋白質與 DAN 中特定序列.

(18) 的結合能力 1、溶劑化自由能 (solvation free energies)17, 20, 25 的預測，. Fig. 2.3 熱力學循環(S. Gentile et al, 2011)7 Fig. 2.3 中，P + L1 代表 P (protein) 與 L1 (ligand 1) 尚未結合的狀態，PL1 則為結合的狀態。其中，ΔGbind(L1), ΔGbind(L2) 為藥物分子與蛋白質結合前後的自由能差異，代表藥物分子與蛋白質之間的結合能力，是我們所關切、期望能加以計算的物理量。但熱力學積分僅能計算相似系統間的自由能差異，必須使用間接的方式求得 ΔGbind(L2) - ΔGbind(L1) 即二藥物分子與蛋白質結合自由能的差異或稱之為 ΔΔGbind 此一物理量可用來衡量二藥物分子與蛋白質結合能力的差異。為求得 ΔΔGbind 先以熱力學積分計算出二相似系統之間的自由能差異即 ΔGfree(L1→L2), ΔGbound(L1→L2) ，由於自由能 G 是一狀態函數與路徑無關，因此經由 2 不同路徑從 P + L1 至 PL2 的 ΔG 會相等： ΔGfree(L1→L2) + ΔGbind(L2) = ΔGbind(L1) + ΔGbound(L1→L2) 加以移項後可得： ΔΔGbind =ΔGbind(L2) - ΔGbind(L1) = ΔGbound(L1→L2) -ΔGfree(L1→L2) 若 ΔΔGbind 為負值代表 ΔGbind(L2) < ΔGbind(L1) ，L2 相較於 L1 與蛋白質有較好的結合能力為較佳的抑制劑分子。反之，若 ΔΔGbind 為正值則代表 L2 與蛋白質的結合能力較差。. 8.

(19) 2.3. Three-transformation approach, TTA 本研究在計算上述的 ΔGfree(L1→L2), ΔGbound(L1→L2) 時使用 TTA7，將二系統. 之間的轉換分為三個步驟 TI1, TI2, TI3 ，TI1 將二系統中相異原子的電荷移除， TI2 進行相異原子的凡得瓦轉換，以及二系統中相同原子的電荷變換，TI3 則重新給予 L2 中相異原子電荷。TTA 的計算，使得二系統中，相異的原子分別經歷了電荷移除 (TI1)、凡得瓦轉換 (TI2)，電荷重新給予 (TI3)，三個步驟，避免在單一原子上同時進行凡得瓦與電荷的轉換。Fig. 2.4 中的數字代表各個原子所帶的電荷，以 –H 轉換為 –OH 為例，說明 TTA 各個步驟中所進行的電荷以及原子種類、數目的轉換。. Fig. 2.4 Three-transformation approach, TTA. 9.

(20) 2.4. 分子模擬系統的建構本研究以 RCSB Protein Data Bank10 資料庫中的 PDB code: 3ZRM 結構為基. 礎，建構分子模擬系統。3ZRM 結構為人類 GSK-3β 蛋白與小分子抑制劑 ZRM 結合的複合體 (complex)，本研究使用 UCSF Chimera 進行分子三度空間結構的觀察與分析，. Fig. 2.5 3ZRM 原始 pdb 構型 Fig. 2.5 中原始的 3ZRM 結構中有 2 個相互對稱的 GSK-3β (A, B chain) 與 ZRM 分子，以及多個結晶水、硫酸根 (SO42-)、甘油 (Glycerol) 分子，其中 B chain 的氨基酸序列較為完整，僅保留 B chain 與 ZRM 分子的複合體以及其中 2 個結晶水進行後續分子模擬系統的建構，如 Fig. 2.6。. 10.

(21) Fig. 2.6 ZRM- GSK-3β (B chain)複合體結晶水的保留：文獻中 8，分析 ZRM 分子與周圍環境氫鍵作用力時提到，其中位於結合區 (binding pocket) 深處的結晶水可能與 ZRM 分子有氫鍵的形成如 Fig. 2.7，考量此一因素將二個結晶水保留。. 11.

(22) Fig. 2.7 ZRM 分子與結晶水、GSK3-β 之間形成氫鍵(G. Gentile et al 2011)8 進一步本研究使用 Gaussian 036 以及 AmberTools11 中的 antechamber, tLeap 進行後續模擬所需檔案的製備。Amber 軟體中對於常見的生物分子，如：蛋白質、脂質、醣類、DNA、RNA 等，內建了完整的作用力場參數，本研究中，使用了 Amber11 內建的 gaff, ff03.r1 二個參數組。但對於小分子藥物我們必須使用其他方式給予分子中電荷的分布，本研究利用 Gaussian03 計算小分子藥物的 RESP (Restrained electrostatic potential) 電荷分布. 30. ，使用 Hartree-Fork 搭配基底函數. 6-31G* (HF/6-31G*)。使用 tLeap 為系統加上 5 個 Cl- 以維持整個系統的電中性。加入 TIP3P 水盒子，模擬系統在水溶液中的狀態，設定 GSK3-β 與 ZRM 複合體與水盒子邊界之最短距離為 12Å 。經 tLeap 處理，產生後續計算所需的參數拓樸學檔案 (.prmtop) 以及原子座標檔案(.inpcrd)。 12.

(23) 使用 DS ViewerPro 6.0 軟體，對 ZRM 分子上的官能基進行轉換成為結構類似物(Mutant)分子，計算 Mutant 分子的 RESP 電荷，使用 tLeap 中 setbox 的功能，定義正確的水盒子尺寸資訊，產生後續計算所需檔案。ZRM 與 GSK3-β 複合體系統，原子總數為 67788。另外，將 ZRM 分子單獨提出，加上水盒子，完成 ZRM 與水溶液系統，原子總數為 3800。同樣使用 DS Viewer 與 tLeap 完成 Mutant 分子的水溶液系統。. 熱力學積分 - 分子動力學. 2.5. (thermodynamic. integration-molecular dynamic, TI-MD) 的模擬計算 2.5.1. 計算步驟. 本研究以 TI-MD 的方式計算二相似系統間自由能的差異，分為 TI1, TI2, TI3 三個步驟，如 Fig. 2.8 各步驟分別進行： 1.. 能量最小化，使用 4000 步最陡下降法 (steepest descent) 以及 1000 步共軛梯度法 (conjugate gradient) ，以找到相對能量較為穩定的構型。其中 TI2 因為凡得瓦轉換的限制，使用 5000 步最陡下降法尋找能量的低點。. 2.. 等體積升溫，在 10 ps 過程中將系統從 10 K 提升至 300 K. 3.. 等壓模擬，分段進行 50 ps, 200 ps, 600 ps (部分計算將 600 ps 延長至 1200 ps) 的模擬計算，收集最後 200 ps 以及 600 ps (or 1200ps)中 dV/dλ 的資訊，計算 ΔG 值。其中，蛋白質與小分子藥物複合體 (complex) 系統中 TI2 的部分，為了使 ΔG 達到良好的收斂，使用了 6 個不同的亂數種子，於計算中對各原子產生不同的初始速度，最後以 6 條軌跡的平均值計算 ΔG。. 13.

(24) Fig. 2.8 TI-MD 計算流程. 2.5.2. 計算相關設定. 等體積升溫的過程中，設定 NTB=1(constant volume periodic boundaries)，並忽略二原子間距離超過 9Å 以上的作用力 (CUT=9)。使用 Langevin dynamics11 碰撞頻率為 2.0ps-1 控制溫度 (NTT=3, gamma_ln=2.0)。等壓模擬過程中，將模擬系統控制在一大氣壓下 (NTB=2)，同樣忽略二原子間距離超過 9Å 以上的作用力 (CUT=9)。在 TI1 與 TI3 的模擬中，各個模擬步驟之間的時間間隔為 2fs (dt=0.002)，為避免氫原子在較大的時間間隔中震盪頻率過大，影響系統的穩定性，使用了 SHAKE19 對氫原子相關的鍵結加入了限制 (NTC=2, NTF=2)。 TI2 的模擬，涉及了二不同系統中原子的消失與生成，本研究使用 Softcore Potential2, 26 定義二不同系統中相異原子的凡得瓦作用位能(ifsc=1)。 TI2 中，因涉及了凡得瓦轉換，無法對二系統中相異進行轉換的原子使用 SHAKE 加以限制氫原子的震盪頻率，在 ZRL → ZRK 以及 ZRM → Mu1~Mu9 的計算中，關閉 SHAKE (NTC=1, NTF=1)，並將時間間隔縮短為 1fs (dt=0.001)。此一設定 TI2 時間間隔為 TI1, TI3 的 1/2，為達到相同的模擬時間 TI2 需進行二倍的模擬步驟，相較於 TI1, TI3 耗費較多的計算資源。而在 ZRM → M15~M24 的計算中，改變計算設定，開啟 SHAKE (NTC=2, NTF=2)，並對二系統中相異進行凡得瓦轉換的原子使用 noshakemask，時間間隔. 14.

(25) 設為 2fs (dt=0.002)。此一設定的結果，使整體系統中除相異原子外，與氫原子相關的鍵結受到了限制，增加系統在較大的時間間隔 2fs 模擬中的穩定性， noshakemask 使二系統中相異原子不受 SHAKE 的限制，得以進行凡得瓦轉換。時間間隔為 2fs 的模擬節省了約 1/2 的計算資源，為更加確認 TI2 中 ΔG 計算結果達到良好的收斂，在 ZRM → M17~M24 中將原先 600ps 的模擬延長為 1200ps。. 15.

(26) Chapter 3 3.1. 結果與討論. 已知抑制劑分子 ZRL, ZRK 間相對結合自由能計算為確認 TI-MD 計算程序的準確程度，首先嘗試計算已知抑制劑分子 ZRL 與. ZRK 之間的相對結合自由能(ΔΔGbind)，與文獻中 8 實驗結果比較。. Fig. 3.1 ZRL → ZRK 熱力學循環 Fig. 3.1 為 ZRL → ZRK 熱力學循環，如 2.2 節中所述，本研究以 TI-MD 的方式計算 ΔGsolution 與 ΔGcomplex 求得 ΔΔGbind，根據 ΔΔGbind 估計二不同小分子藥物 ZRL、 ZRK 與 GSK-3β 結合能力的差異。 ∆∆G𝑏𝑖𝑛𝑑 = ∆G𝑏𝑖𝑛𝑑 (ZRK) − ∆G𝑏𝑖𝑛𝑑 (ZRL) = ∆G𝑐𝑜𝑚𝑝𝑙𝑒𝑥 − ∆G𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 如 2.3 節所述，本研究以 Three-transformation approach, TTA 的方式計算 ΔGsolution 與 ΔGcomplex ，取模擬過程中，最後 600ps 與 200ps 過程中的 dV/dλ 資訊，計算 ΔG 值。在 Table 3.1 中可觀察到，除了 ΔGcomplex TI2 之外，200ps 與 600ps 計算結果差距不大，約在 0.1 kcal/mol 以內。而 ΔGcomplex TI2 則有約 0.3 16.

(27) kcal/mol 的差距，是計算誤差的主要來源。 Table 3.1 ZRL → ZRK 單一軌跡計算結果 (kcal/mol) ΔGcomplex. ΔGsolution. ΔΔGbind. Mutational step 200ps. 600ps. 200ps. 600ps. TI1. 0.65. 0.65. 0.56. 0.57. TI2 (1 traj). 0.27. 0.54. -1.22. -1.12. TI3. 0.92. 0.92. 1.45. 1.44. total. 1.84. 2.12. 0.79. 0.89. 200ps. 600ps. 1.05. 1.23. 為了改善計算誤差，針對 ΔGcomplex TI2 的部分以不同的亂數種子(seeds)，額外計算 5 條軌跡求取平均，結果整理於 Table 3.2。 Table 3.2 ΔGcomplex TI2 不同軌跡計算結果 Seeds. TI2 200ps. TI2 600ps. 71277. 0.27. 0.54. 71379. 0.15. 0.39. 4113. 0.64. 0.51. 275647. 0.40. 0.42. 50111. 0.79. 0.65. 10529. 0.45. 0.49. average. 0.45. 0.50. 為更進一步確認 ΔGcomplex TI2 600ps 模擬過程中的計算誤差，每隔 20ps 計算一次 ΔG 值，對時間作圖如 Fig. 3.2。. 17.

(28) Fig. 3.2 ZRL → ZRK ΔGcomplex TI2 600ps 計算結果與模擬時間關係圖在 Fig. 3.2 中可觀察到，不同軌跡所計算出的 ΔG 值在最初最大達到約 1 kcal/mol 的差距，隨著模擬時間的增加差距逐漸縮減，呈現收斂的現象。深黑色的平均值曲線，漸趨平緩，600ps 模擬過程中的變化約在 0.3 kcal/mol 的範圍之內。 Table 3.3 中以 ΔGcomplex TI2 600ps 6 條軌跡的平均值計算 ΔΔGbind，估計整體的計算誤差約在 0.5 kcal/mol，各計算數值僅取至小數點第一位。. 18.

(29) Table 3.3 ZRL → ZRK ΔΔGbind 600ps 計算結果 600ps. ΔGcomplex. ΔGsolution. TI1. 0.7. 0.6. TI2. 0.5. -1.1. TI3. 0.9. 1.4. total. 2.1. 0.9. ΔΔGbind. 1.2. 文獻中 8 報導，經由實驗測量 ZRL 與 ZRK 的 pIC50 值分別為 6.8 與 6.9，將文獻中的實驗測量值轉換為 ΔΔG ，與計算結果比較。ZRL 與 ZRK 的 pIC50 值，換算為 IC50 值分別為 0.16 與 0.13 μM。本研究中，小分子藥物與 ATP 競爭性的抑制 GSK-3β 活性，可根據 Cheng-Prusoff equation4,. 16. 將 IC50 轉換為 Ki. (inhibition constant, 抑制常數)， 𝐼𝐶50 [𝑆] 1+𝐾 𝑚. 𝐾𝑖 =. 更進一步，根據 Van 't Hoff equation 自由能與 Ki 的關係式， ∆𝐺 = −𝑅𝑇 𝑙𝑛. 1 𝐾𝑖. 將二式應用於 ZRL 與 ZRK ΔΔGbind 的計算， ∆∆𝐺𝑏𝑖𝑛𝑑(𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡) = ∆𝐺𝑏𝑖𝑛𝑑(𝑍𝑅𝐾) − ∆𝐺𝑏𝑖𝑛𝑑(𝑍𝑅𝐿) = −𝑅𝑇 𝑙𝑛. 1 𝐾𝑖(𝑍𝑅𝐾). − (−𝑅𝑇 𝑙𝑛. 1 𝐾𝑖(𝑍𝑅𝐿. = 𝑅𝑇 𝑙𝑛 𝐾𝑖(𝑍𝑅𝐾) − 𝑅𝑇 𝑙𝑛 𝐾𝑖(𝑍𝑅𝐿) = 𝑅𝑇 𝑙𝑛. 𝐾𝑖(𝑍𝑅𝐾) 𝐾𝑖(𝑍𝑅𝐿). 19. ).

(30) 𝐼𝐶50(𝑍𝑅𝐾) [𝑆] 1+𝐾 𝑚 = 𝑅𝑇 𝑙𝑛 𝐼𝐶50(𝑍𝑅𝐿) [𝑆] 1+𝐾 𝑚 = 𝑅𝑇 𝑙𝑛. = 𝑅𝑇 𝑙𝑛. 𝐼𝐶50(𝑍𝑅𝐾) 𝐼𝐶50(𝑍𝑅𝐿). 0.13 0.16. 氣體常數 R = 1.9858775 × 10−3 kcal K−1 mol−1 ，絕對溫度 T = 300K ∆∆𝐺𝑏𝑖𝑛𝑑(𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡) = −0.12 𝑘𝑐𝑎𝑙/𝑚𝑜𝑙 實驗值與 TI-MD 所計算出的 ΔΔGbind 之間約有 1.3 kcal/mol 的差距與近期文獻中 18, 31 以相似 TI-MD 計算流程所得的結果相近。 ZRL → ZRK 的例子中，實驗值與 TI-MD 計算的誤差約為 1.3 kcal/mol， TI-MD 計算誤差約為 0.5 kcal/mol，在後續 3.3 的計算中將嘗試預測未知小分子藥物，本研究中稱之為 Mutants 與 ZRM 之間 ΔGbind 的差異，即 ΔΔGbind(ZRM Mutants)，考量以上二誤差來源，期望能尋找到. →. ΔΔGbind(ZRM → Mutants) 計算結果小於 -2. kcal/mol 的小分子藥物，能較為肯定的推測未知小分子藥物 (Mutants) 相較於 ZRM 對於 GSK-3β 有較佳的抑制能力。. 3.2. ZRL → ZRK 凡得瓦參數調整 Amber11 中對於原子間的凡得瓦作用力以 Lennard-Jones potential 加以描述，. 其中有關各種不同的原子凡得瓦作用力參數多由量子力學的計算而得，嘗試修改 ZRL 中 Br 原子的凡得瓦參數，探討對計算結果的影響。將原始的 ZRL 參數拓樸學檔案 (.prmtop) 中，有關 Br 原子凡得瓦作用力的參數 (LENNARD_JONES_ACOEF, LENNARD_JONES_BCOEF)，調整為原先的 20.

(31) 0.8 倍，以相同的程序進行 TI-MD 的計算。同樣的，為確認、減少計算誤差，在 ΔGcomplex TI2 600ps 模擬過程中，每隔 20ps 計算一次 ΔG 值，對時間作圖如 Fig. 3.3。可觀察到，相較於參數調整前的計算，6 條軌跡的曲線較為分散，但平均值曲線在 300ps 後趨於平緩呈現收斂，可確認計算誤差控制在合理的範圍之內。. Fig. 3.3 ZRL → ZRK ΔGcomplex TI2 600ps 調整凡得瓦參數，計算結果與模擬時間關係圖. 21.

(32) 計算調整凡得瓦參數後的 ΔΔGbind 結果整理於 Table 3.4，在 TI1~TI3 的計算中，TI3 重新給予 ZRK 中變動的原子電荷，與 ZRL 上的 Br 原子無關，此部分沿用先前計算的數值。TI1 的部分，調整參數僅造成微小的改變，ΔΔGbind 的變化主要來自 TI2 由 ZRL 至 ZRK 的凡得瓦轉換。 Table 3.4 ZRL → ZRK ΔΔGbind 600ps 調整凡得瓦參數，計算結果 600ps. ΔGcomplex. ΔGsolution. TI1. 0.7. 0.6. TI2. -0.6. -1.9. TI3. 0.9. 1.4. total. 1.0. 0.1. ΔΔGbind. 0.8. 調整參數後 ΔΔGbind 計算結果為 0.8 kcal/mol ，相較調整前的 1.2 kcal/mol 更為接近 ΔΔGbind(experiment) ，將誤差由原先的 1.3 kcal/mol 縮小至 0.9 kcal/mol。 ZRL → ZRK 計算結果，改變單一原子 Br 的凡得瓦參數對於最後 ΔΔGbind 的計算結果產生了影響，減少了與實驗之間的誤差達 0.4 kcal/mol。若能收集更多蛋白質與小分子藥物結合能力的資訊，藉以修正凡得瓦參數，也許是提高 TI-MD 計算準確性可行的方向。. 22.

(33) 3.3. ZRM, 結構類似物(Mutants)間相對結合自由能計算本節中的計算以已知抑制劑分子 ZRM 為基礎，轉換部分官能基形成 ZRM 的. 結構類似物，本研究中稱之為 Mutants。計算 ZRM 分子與 Mutants 之間對於 GSK-3 的相對自由能，藉以預測 Mutants 相對於 ZRM 之抑制能力。為方便描述進行官能基代換的位置，Fig. 3.4 中根據 IUPAC 命名規則標示 ZRM 分子環上各個原子的編號。ZRM 分子中主結構為圖中左上角，含有 S, N, C=O 的五元、六元融合環(fused ring)部分，為噻吩並吡啶酮(thienopyridinone)，本研究中將此部分簡稱為主結構。主結構 2 號碳位置連接吡啶基(pyridyl group)，7 號碳位置連接羥基苯基(hydroxyphenyl group)。. Fig. 3.4 ZRM 分子結構與原子編號使用 3.1 節 ZRL → ZRK 計算相似的程序，進行 ZRM → Mutants 的計算，本研究將 Mutants 分子命名為 Mu1~M24。Mu1~Mu7 在 ZRM 分子主結構以及吡啶 23.

(34) 基上進行官能基的轉變。將原先的 –H 變更為極性 –OH 與非極性 –CH3, – CH2CH3 官能基，探討此一轉變對於相對結合自由能的影響。Mu9~M24 在 ZRM 分子羥基苯基結構上進行官能基的改變。為更進一步確認、減少 ΔGcomplex TI2 的計算誤差，部分 Mutants 將原先 600ps 的計算時間延長至 1200ps 。在繪製 ΔGcomplex TI2 計算結果與模擬時間關係圖時，為保持一致性，統一將 X 軸的範圍設為 0~1200ps，Y 軸最大值與最小值的差距，除非各軌跡的差異超過 3 kcal/mol ，否則皆訂為 3 kcal/mol ，方便比較不同計算間的差異性。ΔGcomplex TI2 計算結果與模擬時間關係圖以及 ΔΔGbind 600ps 計算結果表格取部分具代表性的結果節錄於本文獻中。. 3.3.1. Mu1. Fig. 3.5 Mu1 分子結構與周遭 GSK-3β 蛋白質環境 Mu1 將 ZRM 主結構中 3 號碳位置的 –H 取代為 –CH3，在 Fig. 3.5 的右圖中以 UCSF Chimera 觀察周遭的蛋白質環境，橘紅色調的區域代表有較高的疏水性，藍色色調的部分代表有較高的親水性。Mu1 於疏水性的區域加入疏水性的 –CH3 官能基，測試其對於 ΔΔGbind 計算結果的影響。 ΔGcomplex TI2 進行了 600ps 6 條軌跡的計算 Fig. 3.6，average 曲線自 200ps 後 24.

(35) 持續上升，相較於先前的計算可能有較大的計算誤差，但估計仍在 0.5 kcal/mol 範圍之內。. Fig. 3.6 ZRM → Mu1 ΔGcomplex TI2 600ps 計算結果與模擬時間關係圖 ZRM → Mu1 ΔΔGbind 計算結果為 0.4 kcal/mol 顯示 Mu1 的抑制能力略低於 ZRM，Table 3.5。 Table 3.5 ZRM → Mu1 ΔΔGbind 600ps 計算結果 600ps. ΔGcomplex. ΔGsolution. TI1. -2.7. -3.1. TI2. -4.2. -3.5. TI3. 1.1. 0.3. total. -5.8. -6.2. 25. ΔΔGbind. 0.4.

(36) 3.3.2. Mu2. Fig. 3.7 Mu1 分子結構與周遭 GSK-3β 蛋白質環境相較於 Mu1，Mu2 在相同位置改以 –OH 取代，ZRM → Mu2 ΔΔGbind 計算結果為 -0.3 kcal/mol 顯示 Mu2 的抑制能力略高於 ZRM。. 3.3.3. Mu3. Fig. 3.8 Mu3 分子結構與周遭 GSK-3β 蛋白質環境 Mu3 於 ZRM 主結構的 6 號碳上以 –CH3 取代原先的 –H ，透過 UCSF Chimera 觀察周遭的蛋白質環境，Mu3 的 –CH3 官能基與 GSK-3β 上 PHE 67 的苯環結構鄰近，同為疏水性的官能基。 ZRM → Mu3 ΔΔGbind 計算結果為 0.2 kcal/mol 顯示 Mu3 的抑制能力略低於 ZRM。 26.

(37) 3.3.4. Mu4. Fig. 3.9 Mu3 分子結構與周遭 GSK-3β 蛋白質環境相較於 Mu3， Mu4 在相同位置改以 –OH 代換，Fig. 3.9 中可觀察到代換上的 –OH 周遭為較疏水的環境，主要由上節中所提到的 PHE 67 的苯環結構所造成。 ZRM → Mu3 ΔΔGbind 計算結果為 1.3 kcal/mol 是目前計算中較為明顯的差異，顯示 Mu4 抑制能力低於 ZRM ，在疏水的周圍環境中代換上親水性的官能基 – OH 可能是造成 Mu4 ΔGbind 上升的主因。. 3.3.5. Mu5. Fig. 3.10 Mu5 分子結構與周遭 GSK-3β 蛋白質環境. 27.

(38) Mu5 於 ZRM 分子吡啶基 (pyridyl group) 結構中 3 號碳位置進行 –OH 的代換。Fig. 3.10 中可觀察到 N 原子與 –OH 分別有機會與鄰近的 GSK-3β 蛋白質骨架 (backbone) 形成氫鍵，或許有助於穩定蛋白質-抑制劑複合體 (protein-inhibiter complex) 的結構。得到較低的結合自由能 (ΔGbind)。計算結果與預期相反，ZRM → Mu5 ΔΔGbind 為 2.2 kcal/mol 顯示 Mu5 之抑制能力明顯低於 ZRM。在這個例子中，可了解到蛋白質-抑制劑複合體中的交互作用十分複雜，以單一固定構型分析推斷結合能力，很可能產生錯誤的偏見。. 3.3.6. Mu7. Fig. 3.11 Mu7 分子結構與周遭 GSK-3β 蛋白質環境 Mu7 分子將原先 Mu1 中 –CH3 的部分改以 –CH2CH3 代換，嘗試在相同位置放入更大的官能基，了解其對於 ΔΔGbind 的影響。. 28.

(39) Fig. 3.12 ZRM → Mu7 ΔGcomplex TI2 600ps 計算結果與模擬時間關係圖 Fig. 3.12 中可觀察到 Mu7 ΔGcomplex TI2 600ps 的計算結果，不同的軌跡彼此歧異、分散的程度略較其他的 Mutants 要來的高。推測 –CH2CH3 官能基深入 GSK-3β 的結合區 (binding pocket) 中，可能與周遭的原子形成空間障礙，使得 MD 模擬的過程中較為不易尋找到穩定的能量低點，造成 ΔG 計算結果隨模擬時間震盪的幅度較高。 Table 3.6 ZRM → Mu7 ΔΔGbind 600ps 計算結果 600ps. ΔGcomplex. ΔGsolution. TI1. -2.7. -3.1. TI2. 2.9. 1.5. TI3. 0.8. 0.3. total. 1.0. -1.3 29. ΔΔGbind. 2.3.

(40) Table 3.6 中 ZRM → Mu7 ΔΔGbind 600ps 計算結果為 2.3 kcal/mol 顯示 Mu7 之抑制能力明顯低於 ZRM ，如前段所述可能由於空間障礙，升高 Mu7 的結合自由能。. 3.3.7. Mu9. Fig. 3.13 Mu9 分子結構與周遭 GSK-3β 蛋白質環境在 Mu9 中，首次嘗試同時於 2 個位置進行官能基的代換。保留 Mu3 中所代換的–CH3，於 ZRM 分子羥基苯基(hydroxyphenyl group)中 3 號碳的位置代換上– OH 官能基，試圖探討多個位置的官能基代換對 ΔΔGbind 計算結果的影響。. 30.

(41) Fig. 3.14 ZRM → Mu9 ΔGcomplex TI2 600ps 計算結果與模擬時間關係圖在 Fig. 3.14 ZRM → Mu9 ΔGcomplex TI2 600ps 的計算中，單以 6 條軌跡計算的結果收斂的情形不甚理想，為確保計算誤差在預期的 0.5 kcal/mol 的範圍內，額外增加了 2 條軌跡的計算。可觀察到與模擬之初各軌跡計算出的 ΔG 值最大差距約為 4 kcal/mol 遠高於先前 Mu1~Mu7 的計算結果。推測，TI-MD 計算中較多原子的轉變可能造成各個軌跡之間較大的歧異，後續 Mu14~Mu24 的計算同樣包含多個官能基的轉換，也有相似的現象，進一步確定了此一推測。. 31.

(42) Table 3.7 ZRM → Mu9 ΔΔGbind 600ps 計算結果 600ps. ΔGcomplex. ΔGsolution. TI1. -12.1. -11.9. TI2. 14.1. 14.1. TI3. -2.7. -2.3. total. -0.7. -0.2. ΔΔGbind. -0.5. Table 3.7 中預測 ΔΔGbind 為 -0.5 kcal/mol 顯示相較於 ZRM， Mu9 有較高的抑制能力，雖然仍未小於期望的 -2.0 kcal/mol ，但相較於 Mu3 羥基苯基 (hydroxyphenyl group)上–OH 的代換使 ΔΔGbind 降低了 0.7 kcal/mol。在後續 M14~M24 的計算中持續探討羥基苯基上不同位向、不同官能基代換所造成的影響。 Table 3.8, Table 3.9 中整理了目前進行的 7 組計算中的分子構型以及計算結果。Mu1~Mu7 中於 ZRM 分子主結構以及吡啶基(pyridyl group)部分進行官能基的轉變，根據 ΔΔGbind 計算結果預測 Mu1~Mu7 分子之抑制能力與 ZRM 沒有明顯差距(Mu1~Mu3)或是明顯低於 ZRM 分子(Mu4~Mu7)。 7 組計算中，我們嘗試根據官能基與周遭環境的疏水、親水、氫鍵作用力推測 ΔΔGbind 計算的結果，但大多未如預期。小分子抑制劑與 GSK-3β 之間的交互作用十分複雜，在模擬過程中彼此的構型交互影響改變，僅以單一構型進行作用力的分析，甚難合理的推測官能基轉變對於抑制能力的影響。 Fig. 3.15 中可以觀察到 ZRM 與 GSK-3β 之結合模式(binding mode)，相較於羥基苯基，ZRM 分子主結構與吡啶基部分位於結合位(binding pocket)較深處，與周遭 GSK-3β 的氨基酸殘跡較為鄰近，在此一部分進行官能基的轉變，可能由於空間 32.

(43) 障礙造成 Mutant 分子抑制能力下降。羥基苯基結構周圍則相對有較多的空間，後續的計算 M14~M24 嘗試於此處進行官能基轉變，期待找到抑制能力明顯優於 ZRM 的 Mutant 分子。. Fig. 3.15 ZRM 分子與 GSK-3β 之結合模式. 33.

(44) Table 3.8 ZRM → Mutants 計算結果 1 Ligand scaffold. Name. Mutants. Gbind (kcal/mol). Mu1. 0.4. Mu2. -0.3. Mu3. 0.2. Mu4. 1.3. ZRM. 34.

(45) Table 3.9 ZRM → Mutants 計算結果 2 Ligand scaffold. Name. Mutants. Gbind (kcal/mol). Mu5. 2.2. Mu7. 2.3. Mu9. -0.5. ZRM. 35.

(46) 3.3.8. ZRM → M14~M24. M14~M24 於 ZRM 分子羥基苯基(hydroxyphenyl group)結構部分進行官能基的轉變，各組計算的 ΔGcomplex TI2 計算結果與模擬時間關係圖，以及 ΔΔGbind 部分具代表性的計算結果整理於 Fig. 3.18~Fig. 3.21 與 Table 3.11~Table 3.14 中， M14~M24 包含多個官能基的轉變，不同軌跡 ΔGcomplex TI2 計算結果有較大的差異，但 average 曲線在 600ps 或 1200ps 的計算中皆呈現合理收斂，可以確認計算誤差在預期範圍之內。Table 3.10 中整理 M14~M24 分子構型，以及 ΔΔGbind 計算結果。探討 Mu3, Mu9 與 M24 計算結果的比較，ΔΔGbind 分別為 0.2, -0.5, 0.5 kcal/mol。Fig. 3.16 為 3 組 Mutant 分子結構的比較。Mu3 在原 ZRM 分子主結構 6 號碳的位置增加了 –CH3 官能基，M21 則在原 ZRM 分子羥基苯基結構中原本已存在 –OH 的鄰位多增加了一個 –OH 官能基，Mu9 則兼具上述二官能基的轉變。單獨增加 –CH3, –OH 官能基 ΔΔGbind 分別為 0.2, 0.5 kcal/mol ，顯示單獨增加二官能基對抑制能力有負面的影響。而 Mu9 中同時增加二官能基，ΔΔGbind 則達到 -0.5 kcal/mol 顯示同時增加二官能基對抑制能力有正面的影響。綜合以上結果，可了解到此一例子中個別官能基的負面影響不具有加成性，小分子藥物中官能基的轉變與抑制能力的關係不易以個別的計算結果加以預測。. Fig. 3.16 Mu3, Mu9 與 M24 分子結構 36.

(47) Table 3.10 中 M19, M21 計算結果有較為突出的表現，ΔΔGbind 分別為 -0.6, -0.8 kcal/mol 顯示此二 Mutant 分子相較於 ZRM 對於 GSK-3β 有較高的抑制能力。. Fig. 3.17 M19 與 M21 分子結構 Fig. 3.17 中可觀察到 M19 與 M21 分子結構十分相近，M21 相較 M19 多了一個 –F 官能基，二分子同樣在原 ZRM 分子的羥基苯基結構中 2, 5 號碳的位置增加了二個互為對位的 –OH 官能基。此特定位置官能基的轉變，似乎有助於提升 Mutant 分子對於 GSK-3β 的抑制能力。使用 AmberTool11 中的 ptraj 程式，分析計算結果最佳的 M21 分子於 1200ps 模擬過程中的軌跡。. 37.

(48) Table 3.10 Ligand scaffold. ZRM → Mutants 計算結果 3. Name. R groups. Gbind (kcal/mol). M14. 1.0. M15. 1.3. M17. 0.9. M19. -0.6. M21. -0.8. M23. 0.5. M24. 0.5. ZRM. 38.

(49) Fig. 3.18 ZRM → M17 ΔGcomplex TI2 1200ps 計算結果與模擬時間關係圖. Table 3.11. ZRM → M17 ΔΔGbind 1200ps 計算結果. 1200ps. ΔGcomplex. ΔGsolution. TI1. -7.1. -6.8. TI2. 29.5. 28.3. TI3. -4.1. -3.8. total. 18.3. 17.6. 39. ΔΔGbind. 0.7.

(50) Fig. 3.19 ZRM → M19 ΔGcomplex TI2 1200ps 計算結果與模擬時間關係圖. Table 3.12. ZRM → M19 ΔΔGbind 1200ps 計算結果. 1200ps. ΔGcomplex. ΔGsolution. TI1. 20.4. 20.3. TI2. 5.0. 4.8. TI3. -42.0. -41.1. total. -16.6. -16.0. 40. ΔΔGbind. -0.6.

(51) Fig. 3.20 ZRM → M21 ΔGcomplex TI2 1200ps 計算結果與模擬時間關係圖. Table 3.13. ZRM → M21 ΔΔGbind 1200ps 計算結果. 1200ps. ΔGcomplex. ΔGsolution. TI1. 20.4. 20.3. TI2. 5.0. 4.8. TI3. -42.0. -41.1. total. -16.6. -16.0. 41. ΔΔGbind. -0.6.

(52) Fig. 3.21 ZRM → M23 ΔGcomplex TI2 1200ps 計算結果與模擬時間關係圖. Table 3.14. ZRM → M23 ΔΔGbind 1200ps 計算結果. 1200ps. ΔGcomplex. ΔGsolution. TI1. -7.1. -6.8. TI2. 30.3. 29.7. TI3. -11.7. -11.8. total. 11.5. 11.1. 42. ΔΔGbind. 0.5.

(53) ZRM, M21 與 GSK-3β 複合結構於 1200ps 模擬過程. 3.4. 中的軌跡分析。本節中對 GSK-3β- M21 複合結構進行 MD 模擬。模擬步驟如 2.5.1 節中所述包含了：1.能量最小化 2.等體積升溫 3.等壓模擬三個部分，分析最後 1200ps 過程中系統原子的軌跡。. 3.4.1. 羥基苯基與主結構之間二面角 (dihedral angle) 分析. 在本節中，期望探討羥基苯基在模擬過程中是否會有翻轉的現象。如 Fig. 3.22 所示，羥基苯基與主結構之間以單鍵鍵結，在一般情況下應可自由旋轉，A, B 二構型可視為相同的分子。但在 M21 與 GSK-3β 複合結構中，抑制劑分子與蛋白質之間的交互作用可能限制了此一部份結構的轉動，A, B 二構型的轉換必須跨越一定的位能障礙 (potential barrier)。在本研究 200ps 與 1200ps 模擬過程中 M21 分子能否越過位能障礙可透過 AmberTool11 中的 ptraj 程式加以確認。. A. B. Fig. 3.22 M21 分子結構 Fig. 3.23 為 M21 分子結構與二面角示意圖，使用 ptraj 程式分析 1200ps 模擬過程中，C1, C11, C17, C16 四個原子所形成的二面角 (dihedral angle) 在模擬過程中的變化，以同樣方式分析 ZRM 分子進行比較。以二面角對時間作圖如 Fig. 43.

(54) 3.24。可觀察到 M21 的結構相對穩定，二面角約在 40° 上下震盪。而 ZRM 則似乎有二個相對的能量低點，二面角約為 30° 與 -30°，分子構型在二者之間震盪、轉換。二面角分析的結果顯示，M21 分子中的羥基苯基結構在 1200ps 模擬時間中，無法進行翻轉，Fig. 3.22 中 A, B 二構型的轉變無法在本研究所使用的計算程序中達成。因此，在此一計算程序下 A, B 二構型應視為不相同的二個分子，預期二者的計算結果也會有一定程度的差異。. Fig. 3.23 M21 分子結構與二面角示意圖. 44.

(55) Fig. 3.24 ZRM, M21 與 GSK-3β 複合結構中，二面角與時間關係圖。. 3.4.2. M21 分子與周圍原子的氫鍵分析. 使用 AmberTool11 中的 ptraj 程式，進行氫鍵分析。在氫鍵分析中，設定形成氫鍵的三個原子，如 O—H ᠁O 中三原子所形成的角度必須大於 120°，二個 O 原子之間的距離必須小於 3.5 Å 。 Fig. 3.25 中藍色虛線代表分子間氫鍵，並標示出分子或殘基的名稱，以及在 1200ps 模擬過程中符合氫鍵形成條件的比例。M21 分子與 2.4 節中所保留的水分子 WAT 356 形成穩定的氫鍵，53.00% 與 24.67% 分別代表水分子中二組 O—H 與 M21 上的雙鍵 O 形成氫鍵的百分比。WAT 356 分子與 M21 形成穩定的氫鍵與文獻中 8 以結晶構型預測的結果相符。. 45.

(56) Fig. 3.25 M21 於 1200ps 模擬中與周圍原子形成氫鍵的比例更進一步，在 Fig. 3.26 中針對保留於 GSK-3β 結合位深處的二水分子與周圍原子於 1200ps 模擬中的氫鍵形成百分比進行分析。可觀察到二水分子之間與周遭的 ASP 165, GLU 62 以及 M21 在模擬過程中有極高的比例形成氫鍵，可能有助於穩定此部份構型。結晶水的保留應能提升模擬計算的正確性。在 Fig. 3.25 中，M21 與周遭的氨基酸殘基形成氫鍵，吡啶基中的 N 原子與 VAL 100 形成穩定的氫鍵 (99.83%)，與 Fig. 2.78 中預測的結果相符 (VAL 100 即 Fig. 2.7 中的 V135，分子模型經 tLeap 處理後胺基酸重新編號造成此一差異)。ASP 165 中的二個 O 原子交替形成氫鍵 (70.42%, 27.25%) 則是文獻中所未報導的部分。. 46.

(57) Fig. 3.26 蛋白質結合位深處水分子與周遭原子氫鍵形成的比例值得注意的是 M21 分子中羥基苯基上的 –OH 官能基，與 TIP3P 加入的水分子 WAT 6194 形成穩定的氫鍵 (97.67%)。Fig. 3.27 中以 WAT 6194 為中心，分析其與周遭原子的氫鍵形成比例。WAT 6194 類似橋梁，分別與 M21 以及 GLN 150 之間形成穩定的氫鍵，穩定三者之間的結構。M21 中的 –OH 也許藉此方式提升 M21 與 GSK-3β 複合結構的穩定度，進而提高分子的抑制能力。. 47.

(58) Fig. 3.27 WAT 6194 與 M21 以及 GLN 150 之間氫鍵形成的比例. M19 與 GSK-3β 複合結構構型分析. 3.5. M19 與 M21 構型相似，ΔΔGbind 計算結果分別為 -0.6, -0.8 kcal/mol ，顯示二分子對於 GSK-3β 有較佳的抑制能力，分析 TI-MD 模擬中 200ps, 1200ps 計算過程最後一步驟所儲存的軌跡檔案，檢視 M19 是否具有與 M21 相似的 Mutant-水分子 -GLN 150 穩定架構。 TI-MD 計算中，TI3 λ=0.9 的分子動力學模擬最為接近 M19 與 GSK-3β 複合結構的系統，使用 Amber11 中 ambpdb 指令將上述模擬中 .inpcrd .prmtop 檔案轉為 .pdb 檔案，以 UCSF Chimera 中 FindHBond 功能尋找 M19 與周圍分子的氫鍵。 Fig. 3.28 為 1200ps 計算中最後一步構型分析的結果，與 M21 相似 WAT 5091 分別與 GLN 150 以及 M19 形成氫鍵。Fig. 3.29 則為 200ps 計算結果的分析，M19 直接與 GLN 150 形成氫鍵，WAT 5091 在此一階段位於較遠的位置。 48.

(59) 推測 WAT 5091 在 1200ps 模擬過程中，從較遠處進入 M19 與 GLN 150 之間，使原先 200ps 模擬結果中 GLN 150 與 M19 之間直接的氫鍵鍵結，轉變為 M19-水分子-GLN 150 結構。M19 氫鍵結構的分析，更進一步說明了 Mutant-水分子-GLN 150 結構可能是提升抑制能力的重要因子。. Fig. 3.28 M19 TI3 λ=0.9 1200ps 計算結果氫鍵分析. 49.

(60) Fig. 3.29 M19 TI3 λ=0.9 200ps 計算結果氫鍵分析. 50.

(61) ZRM, 結構類似物(Mutants)間相對結合自由能計算結. 3.6. 果整理本節中，根據 3.3 計算結果將 Mutants 依照抑制能力分類，歸納分子構型與抑制能力之間的關係。. 3.6.1. 抑制能力明顯低於 ZRM，ΔΔGbind > 0.5 kcal/mol. Table 3.15 中整理了 ΔΔGbind > 0.5 kcal/mol 的 Mutant 分子，Mu4, Mu5, Mu7 於 ZRM 分子主結構以及吡啶基部分進行官能基轉換，推測此部分分子周圍空間障礙可能是造成抑制能力下降的原因之一。M14, M15 則於羥基苯基部分進行官能基轉換，相較於 M19, M21 於不同的位向進行 –OH 的官能基轉換，對抑制能力可能有相反的影響。. 3.6.2. 抑制能力略低於 ZRM，0 < ΔΔGbind ≦0.5 kcal/mol. Table 3.16 中的 Mutant 分子 ΔΔGbind 計算結果介於 0~0.5 kcal/mol 之間，於結合位深處或是羥基苯基部分進行官能基變換的分子皆有部分落於此一區間。. 3.6.3. 抑制能力略高於 ZRM，-0.5 ≦ ΔΔGbind < 0 kcal/mol. Table 3.17 中的 Mutant 分子 ΔΔGbind 計算結果介於 0~0.5 kcal/mol 之間。其中，Mu9 嘗試同時於主結構以及羥基苯基部分進行官能基變換，計算結果位於 -0.5 kcal/mol。. 3.6.4. 抑制能力明顯高於 ZRM，ΔΔGbind ≦ -0.5 kcal/mol. Table 3.18 中的 Mutant 分子 ΔΔGbind 計算結果小於 -0.5 kcal/mol。在 3.4 以及 3.5 節中 M19, M21 MD 模擬軌跡以及構型分析有相似的結果，羥基苯基上 2,5 號碳位置的官能基轉變，可能是影響抑制能力的重要因子。 51.

(62) Table 3.15 Ligand scaffold. ΔΔGbind > 0.5 kcal/mol 之 Mutant 分子 Name. Mutants. Gbind (kcal/mol). Mu4. 1.3. Mu5. 2.2. Mu7. 2.3. ZRM. R groups. M14. 1.0. M15. 1.3. M17. 0.9. ZRM. 52.

(63) Table 3.16 Ligand scaffold. 0 < ΔΔGbind ≦0.5 kcal/mol 之 Mutant 分子 Name. Mutants. Gbind (kcal/mol). Mu1. 0.4. Mu3. 0.2. ZRM. R groups. M23. 0.5. M24. 0.5. ZRM. 53.

(64) Table 3.17 Ligand scaffold. -0.5 ≦ ΔΔGbind < 0 kcal/mol 之 Mutant 分子 Name. Gbind. Mutants. (kcal/mol). Mu2. -0.3. Mu9. -0.5. ZRM. Table 3.18 Ligand scaffold. ΔΔGbind < -0.5 kcal/mol 之 Mutant 分子 Name. R groups. Gbind (kcal/mol). M19. -0.6. M21. -0.8. ZRM. 54.

(65) 3.7. 22A, 結構類似物(Mutants)間相對結合自由能計算 22A 為文獻中 8 報導的已知抑制劑分子 pIC50 為 8.0，抑制能力優於 ZRM 的 6.9。. 以 3.1 中的式子計算 ΔΔGbind(ZRM → 22A,experiment)約為 -1.5 kcal/mol。以 22A 分子結構為基礎，在原先含 N 以及 –Cl 官能基的六元雜環上增加 1~2 個 –OH 官能基。以 2.5 中的程序計算 ΔΔGbind(22A →. Mutants). 期望找到抑制能力優於 22A 的分子。各. 組計算的 ΔGcomplex TI2 計算結果與模擬時間關係圖，以及 ΔΔGbind 計算結果整理於附錄中。 Table 3.19 中整理了 Mutants 分子的構型以及計算結果。M11 與 M12 Gbind 計算結果分別為 0.3, -0.1 kcal/mol ，Mutants 與 22A 對於 GSK-3β 之抑制能力並無顯著的差異。 Table 3.19 Ligand scaffold. 22A → Mutants 計算結果. Name. R group. Gbind (kcal/mol). M11. 0.3. M12. -0.1. 22A. 55.

(66) Chapter 4. 結論. 本研究使用熱力學積分分子動力學計算 18 組類似的小分子抑制劑與 GSK-3β 激酶的相對結合自由能。首先，在 ZRL → ZRK 的計算之中，了解到以本研究所使用之 TI-MD 程序計算小分子藥物與 GSK-3β 之相對結合自由能，計算誤差不確定度(error bar)估計為 0.5 kcal/mol。與實驗結果比較有 1.3 kcal/mol 的誤差。在預測結構類似物(Mutants)的結合能上，M21 是本研究計算中相較 ZRM 與 GSK-3β 結合能力最佳的分子 ΔΔGbind 為 -0.8 kcal/mol ，M21 之分子架構可作為後續模擬計算的重要參考。分析 M21 與 GSK-3β 複合體 MD 模擬的軌跡，了解到分子中的羥基苯基結構在 1200ps 模擬中並無翻轉的現象，翻轉前後之構型在本研究之計算程序中應視為不同的分子。分析模擬過程中分子間氫鍵的形成，大致與先前文獻 8 中以單一結晶構型預測的結果一致。於結合位深處保留的二個結晶水分子和周圍環境有穩定的氫鍵形成，結晶水的保留應有助於提升模擬的正確性。 M21 的計算中發現，有一以水分子為橋梁的氫鍵 -由 TIP3P 水分子與 M21 分子上增加的 –OH 以及 GSK-3β 之胺基酸殘基 GLN 150 形成氫鍵-，這形成 M21水分子-GLN 150 的穩定架構，或許是提升 M21 與 GSK-3β 之結合能力的重要因子。整體來說，在本研究計算結果中 Mutant 分子於 ZRM 分子主結構以及吡啶基位置(結合位深處)進行官能基代換，無法得到結合能力明顯高於 ZRM 的分子，部分 Mutant 分子抑制能力明顯低於 ZRM。結合位深處的官能基轉變可能由於周圍的空間障礙，無法明顯的提升分子對於 GSK-3β 的抑制能力。以單一構型觀察抑制劑. 56.

(67) 分子與周圍蛋白質環境的交互作用，推測官能基代換對於 ΔΔGbind 計算結果的影響，絕大多數無法得到預期的結果。於 ZRM 分子羥基苯基結構上進行官能基的轉換，有較高的機會得到抑制能力高於 ZRM 的分子。在 M19, M21 的例子中於羥基苯基 2,5 號碳上加入 –OH 官能基皆有不錯的計算結果。值得以此分子結構為基礎更深入探討，並透過實驗加以印證。. 57.

(68) REFERENCES 1. F. R. Beierlein, G. G. Kneale, and T. Clark, 'Predicting the Effects of Basepair Mutations in DNA-Protein Complexes by Thermodynamic Integration', Biophys J, 101 (2011), 1130-8.. 2. F. P. Buelens, and H. Grubmuller, 'Linear-Scaling Soft-Core Scheme for Alchemical Free Energy Calculations', J Comput Chem, 33 (2012), 25-33.. 3. D. A. Case, T. Darden, T. E. Cheatham, C. Simmerling, J. Wang, R. E. Duke, R. Luo, R. C. Walker, W. Zhang, K. M. Merz, B. Roberts, B. Wang, S. Hayik, A. Roitberg, G. Seabra, I. Kolossvary, K. F. Wong, F. Paesani, J. Vanicek, X. Wu, S. R. Brozell, T. Steinbrecher, H. Gohlke, Q Cai, X. Ye, J. Wang, M. -J. Hsieh, G. Cui, D. R. Roe, D. H. Mathews, M. G. Seetin, C. Sagui, V. Babin, T. Luchko, S. Vusarov, A. Kovalenko, and P. A. Kollman, Amber 11 (San Francisco: University of California, 2010).. 4. Y. Cheng, and W. H. Prusoff, 'Relationship between the Inhibition Constant (Ki) and the Concentration of Inhibitor Which Causes 50 Per Cent Inhibition (I50) of an Enzymatic Reaction', Biochem Pharmacol, 22 (1973), 3099.. 5. N. Embi, D. RYLATT, and P. COHEN, 'Glycogen Synthase Kinase‐3 from Rabbit Skeletal Muscle', Eur J Biochem, 107 (1980), 519-27.. 6. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, T. Vreven Jr., K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M M. Hada, and K. Toyota Ehara, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, C. Adamo, J. 58.

(69) Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez, and J. A. Pople, 'Gaussian 03, Revision A.1', Gaussian, Inc., 2003). 7. S. Genheden, I. Nilsson, and U. Ryde, 'Binding Affinities of Factor Xa Inhibitors Estimated by Thermodynamic Integration and Mm/Gbsa', J Chem Inf Model, 51 (2011), 947-58.. 8. G. Gentile, G. Bernasconi, A. Pozzan, G. Merlo, P. Marzorati, P. Bamborough, B. Bax, A. Bridges, C. Brough, P. Carter, G. Cutler, M. Neu, and M. Takada, 'Identification of 2-(4-Pyridyl)Thienopyridinones as Gsk-3beta Inhibitors', Bioorg Med Chem Lett, 21 (2011), 4823-7.. 9. L. A. Gil, P. A. Valiente, P. G. Pascutti, and T. Pons, 'Computational Perspectives into Plasmepsins Structure-Function Relationship: Implications to Inhibitors Design', J Trop Med, 2011 (2011), 657483.. 10. J. Westbrook H. M. Berman, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, Ilya N. Shindyalov and P. E. Bourne, 'The Protein Data Bank', Nucleic Acids Res, 28 (2000), 235-42.. 11. J. A. Izaguirre, D. P. Catarello, J. M. Wozniak, and R. D. Skeel, 'Langevin Stabilization of Molecular Dynamics', J Chem Phys, 114 (2001), 2090.. 12. F. Wandosell J. Avila, F. Hernández, 'Role of Glycogen Synthase Kinase-3 in 59.

(70) Alzheimer’s Disease Pathogenesis and Glycogen Synthase Kinase-3 Inhibitors', Expert Rev. Neurother., 10 (2010), 703-10. 13. R. S. Jope, and G. V. Johnson, 'The Glamour and Gloom of Glycogen Synthase Kinase-3', Trends Biochem Sci, 29 (2004), 95-102.. 14. M. W. Karaman, S. Herrgard, D. K. Treiber, P. Gallant, C. E. Atteridge, B. T. Campbell, K. W. Chan, P. Ciceri, M. I. Davis, P. T. Edeen, R. Faraoni, M. Floyd, J. P. Hunt, D. J. Lockhart, Z. V. Milanov, M. J. Morrison, G. Pallares, H. K. Patel, S. Pritchard, L. M. Wodicka, and P. P. Zarrinkar, 'A Quantitative Analysis of Kinase Inhibitor Selectivity', Nat Biotech, 26 (2008), 127-32.. 15. J. G. Kirkwood, 'Statistical Mechanics of Fluid Mixtures', J Chem Phys, 3 (1935), 300.. 16. S. Lazareno, and N.J.M. Birdsall, 'Estimation of Competitive Antagonist Affinity from Functional Inhibition Curves Using the Gaddum, Schild and Cheng-Prusoff Equations', Br J pharmacol, 109 (1993), 1110-19.. 17. S. A. Martins, and S. F. Sousa, 'Comparative Assessment of Computational Methods for the Determination of Solvation Free Energies in Alcohol-Based Molecules', J Comput Chem, 34 (2013), 1354-62.. 18. J. Michel , N. Foloppe, and J. W. Essex, 'Rigorous Free Energy Calculations in Structure-Based Drug Design', Molecular Informatics, 29 (2010), 570-78.. 19. S. Miyamoto, and P. A. Kollman, 'Settle: An Analytical Version of the Shake and Rattle Algorithm for Rigid Water Models', J Comput Chem, 13 (1992), 952-62.. 20. D. L Mobley, C. I. Bayly, M. D. Cooper, M. R. Shirts, and K. A. Dill, 'Small Molecule Hydration Free Energies in Explicit Solvent: An Extensive Test of Fixed-Charge Atomistic Simulations', J Chem Theory Comput, 5 (2009), 350-58.. 21. D. L. Mobley, and K. A. Dill, 'Binding of Small-Molecule Ligands to Proteins: 60.

(71) "What You See" Is Not Always "What You Get"', Structure, 17 (2009), 489-98. 22. 'Pharmaceutical Industry Profile 2012',. (Washington, DC: PhRMA, 2012).. 23. S. Phukan, V. S. Babu, A. Kannoji, R. Hariharan, and V. N. Balaji, 'Gsk3beta: Role in Therapeutic Landscape and Development of Modulators', Br J Pharmacol, 160 (2010), 1-19.. 24. R. Seger, and E. G. Krebs, 'The Mapk Signaling Cascade', FASEB J, 9 (1995), 726-35.. 25. T. Steinbrecher, I. Joung, and D. A. Case, 'Soft-Core Potentials in Thermodynamic Integration: Comparing One- and Two-Step Transformations', J Comput Chem, 32 (2011), 3253-63.. 26. T. Steinbrecher, D. L. Mobley, and D. A. Case, 'Nonlinear Scaling Schemes for Lennard-Jones Interactions in Free Energy Calculations', J Chem Phys, 127 (2007), 214108.. 27. C. Sutherland, 'What Are the Bona Fide Gsk3 Substrates?', Int J Alzheimers Dis, 2011 (2011), 505607.. 28. I. V. Tetko, P. Bruneau, H. W. Mewes, D. C. Rohrer, and G. I. Poda, 'Can We Estimate the Accuracy of Adme-Tox Predictions?', Drug Discov Today, 11 (2006), 700-7.. 29. M. Utton U. Wagner, J. Gallo and C. C. J. Miller, 'Cellular Phosphorylation of Tau by Gsk-3b Influences Tau Binding to Microtubules and Microtubule Organisation', Journal of Cell Science, 109 (1996), 1537-43.. 30. R. J. Woods, and R. Chappelle, 'Restrained Electrostatic Potential Atomic Partial Charges for Condensed-Phase Simulations of Carbohydrates', J Mol Struct, 527 (2000), 149-56.. 31. K. W. Wu, P. C. Chen, J. Wang, and Y. C. Sun, 'Computation of Relative Binding 61.

(72) Free Energy for an Inhibitor and Its Analogs Binding with Erk Kinase Using Thermodynamic Integration Md Simulation', J Comput Aided Mol Des, 26 (2012), 1159-69.. 62.

(73) 附錄. 附圖 1. 附圖 2. 藥物開發時程以及所需消耗經費(PhRMA, 2012)22. 不同自由能計算方法之理論精確度以及所需耗費的計算時間(L. A. Gil et al, 2011)9. 63.

(74)