細胞外調節激酶、纖維母細胞生長受體激酶與其抑制劑結合自由能之計算：熱力學積分分子動態及分子嵌合模擬研究

全文

(1)國立臺灣師範大學化學系碩士論文 Department (Graduate institute) of Chemistry, National Taiwan Normal University Master Thesis. 指導教授：孫英傑博士 Advisor: Ying-Chieh Sun, Ph.D.. 細胞外調節激酶、纖維母細胞生長受體激酶與其抑制劑結合自由能之計算：熱力學積分分子動態及分子嵌合模擬研究. Computational Study of Binding Energy of Protein-Ligand Complexes for Two Kinases: Thermodynamic Integration Molecular Dynamics and Docking Simulation 研究生：吳冠緯 Kuan-Wei Wu 中華民國壹零壹年陸月 June, 2012.

(2) 2.

(3) 誌謝師大六年，要感謝的人太多。回想起來，在這不算短的日子裡，自己從一個平凡的高中生，轉變成一個更趨成熟的人。能完成這篇論文，首先我要感謝我的指導教授孫英傑老師帶領我進入分子模擬的世界。自從大二升大三的暑假開始在實驗室裡做專題生後就受到老師許多在生活上以及科學研究上的啟發。依然記得許多老師您與我分享的寶貴人生經驗以及對化學研究的深度並且常常惦記著當初大學時老師對我的鼓勵、分享做學問的態度，以及待人處事的種種。除此之外，我要感謝口試委員中央研究院原子與分子研究所的楊大衍教授以及陽明大學生命科學系暨基因體研究所的許世宜教授，因為您對這份論文的悉心建議及指教，使的我在這兩年內的研究成果能夠更趨完善。特別感謝國家科學委員會的經費提供及國家高速網路計算中心的化學組的何智雄博士，因為有何博士的協助，幫我們解決了許多難以克服的技術問題。大學四年裡，我要感謝我的大學同學凱傑、智偉、修安、瑩璇、士哲、彥凱、瑾怡、玟淑、禺彤；鋰家國和學長、占杰學長、莙淳學姊、雅凡學姊、學伴哲瑋、學弟妹家榮、元錫、界志、知敏。在大學裡認識你們令我的大學生活不虛此行，從你們的故事中我學到了許多寶貴的生活經驗。研究所兩年，除了感謝實驗室前博士後研究員王軍學長、鼎欣學長、思琪學姊、昭同學長、書吉學長，同學博晉、彥樓、永富、孟璇以及學弟興洲、育劭每天生活的分享、研究上的討論切磋。還要感謝在台大念碩班的大學同學乃華、亭宇；在師大念碩班的延展、家睿、宥廷、展輝、岳樺、耘政。你們的鼓勵和串門子，是我碩士班生活最溫馨的一頁。高中同學珈銘、睿珊、柏彥、于萱以及 Grand Toastermasters、台北歌德學院的大家、化學系的助教群以及曾經與我分享許多想法的老師們與朋友們，謝謝你們讓我的碩士生活有不同的啟發。最後，感謝我親愛的爸爸媽媽以及老弟以及所有的親人，因為你們一路的照顧、栽培以及分享，讓我能夠順利完成這本論文。中華民國一零一年六月，于師大分部公館校區. 1.

(4) 2.

(5) 中文摘要癌症一直為十大死因之一。因此與癌症有關的標靶藥物研發是目前國際間藥物化學家研究的主要方向之一。由於在實驗上合成藥物需要大量的人力、物力。因此若可以使用電腦來輔助藥物設計、藥物探索，不僅可以縮短研發的時程，更可以大量減少有機合成研發的成本。蛋白質激酶目前被認為與許多癌症有關。在本篇論文中，我們使用電腦輔助的方法來探討化學分子與蛋白質激酶間的作用力型態以及如何提升篩選具有生物活性的小分子。我們使用兩種理論方法，分別是熱力學積分及分子嵌合。. 在熱力學積分中，我們所針對的蛋白質標的為 ERK2。我們以原本已知的活性分子開始。預測了一系列新的類似物，這些類似物是以苯環上的極性或非極性官能基做為區別。在我們的結果中，羥基(OH) 取代後發現其與蛋白質間的結合能力約比之前研究中的抑制劑分子好 2 kcal/mol，而其他官能基的取代約增進 1 kcal/mol。我們的預測結果可為實驗學家提供有機合成、蛋白質與細胞實驗的理論依據。. 除了分子動力學之外，我們也以分子嵌合計算探討 FGFR1 蛋白質幾個效應對豐富指數(enrichment factor)的影響。豐富指數為評價虛擬篩選優劣的指標之一。在豐富指數的研究中，我們採用 DUD 資料庫中活性及非活性的分子及設定不同的蛋白質條件。在設定條件時發現擺動 514 號胺基酸(賴氨酸)能夠大幅度增加正確挑選出活性分子的數量。並且設定絞鍊片段上 564 號胺基酸(丙胺酸)的氫鍵限制以進一步探究篩選效率。從本論文的研究中，我們尋找到最佳的條件，也希望這些最佳化條件可以進一步成為此類蛋白激酶大量虛擬篩選中的基礎。進而幫助先導化合物的找尋。. 關鍵字: 分子動力學模擬、熱力學積分、分子嵌合、蛋白質激酶. 3.

(6) 4.

(7) ABSTRACT Cancer has been one of the top ten leading causes of death for several decades. The target drugs research for cancer therapy is now a popular field among international medicinal chemists. Because of the significant amount of money and human resources spent in the drug development process, computer-aided drug design method is an attractive tool to reduce cost and assist drug discovery. Protein kinases are one of the protein families which are drug targets for cancer therapy. Here, we selected two kinases, which are ERK2 and FGFR1 kinases, and used computer modeling to investigate binding energy of inhibitor-protein complexes for these two kinases.. In the part of ERK2, we used thermodynamic integration MD method to compute relative binding free energy of several ERK2-inhibitor complexes of interest. We carried out computations to predict G for new analogs, focusing on placing polar and nonpolar functional groups at the meta site of benzene ring, to see if these ligands have better binding affinity than the above ligands. The computations resulted that a ligand with polar –OH group has better binding affinity than the previous examined ligand by ~2.0 kcal/mol and two other ligands have better affinity by ~1.0 kcal/mol. The predicted better inhibitors of this kind should be of interest to experimentalists for future experimental enzyme and/or cell assays.. In addition to TI-MD simulation, we also worked on interactions of FGFR1 kinase-inhibitor complexes using docking computation, focusing on how enrichment factor (EF) enhances in virtual screening by including side chain movement and applying hydrogen bond constraint for this kinase. To this end, active and decoy compounds from the Directory of Useful Decoys. 1. database was obtained and benchmarked with GOLD. program. Interestingly, among combinations of side chains which were allowed to move, EF is significantly higher with movement of Lys514 compared with others. In addition, the effect of adding hydrogen bond constraint at a residue located in the hinge segment, Ala564, was also examined. The results were analyzed and discussed. The present results should be useful for virtual screening of large databases against this kinase.. Key words: Molecular Dynamics、Thermodynamic Integration、Molecular Docking、 Protein Kinase 5.

(8) 6.

(9) 目錄(CONTENTS) 口試委員會審定書 ....................................................................................................... # 誌謝 .............................................................................................................................. 1 中文摘要 ...................................................................................................................... 3 ABSTRACT ................................................................................................................. 5 目錄(CONTENTS) ...................................................................................................... 7 圖目錄 (LIST OF FIGURES) ................................................................................... 11 表目錄 (LIST OF TABLES) ..................................................................................... 13 符號(SYMBOLS)....................................................................................................... 14 第一部分分子動力學模擬 (Part1: Molecular Dynamics Modeling) ................................................................... 16 Chapter 1. 緒論 (Introduction) ............................................................................ 18. 1.1. 前言 (Preface) ......................................................................................... 18. 1.2. 細胞外調節激酶二與癌症的關聯 (The Relationship Between Extracellular Signal Regulated Kinase Protein Kinase 2 and Cancer) . 18. 1.3 Chapter 2 2.1. 研究目標 (Research Objective) .............................................................. 21 模擬方法 (Simulation Methods) ........................................................ 22 理論 .......................................................................................................... 22. 2.1.1. 分子力學與分子動力學 ......................................................................... 22. 2.1.2. 自由能計算 ............................................................................................ 24. 2.1.3. 熱力學循環與熱力學積分 ..................................................................... 25. 2.2. 模型設定 (Modeling Setup) .................................................................... 27. 2.2.1. 溶液的分子 ............................................................................................ 27. 2.2.2. 蛋白質中的分子 .................................................................................... 27. 2.2.3. 模擬參數 ................................................................................................ 28. Chapter 3. 結果與討論 (Results and Discussions) .............................................. 30. 3.1. CF3 變換為 82A ...................................................................................... 33. 3.2. OHW 變換為 82A.................................................................................... 36. 3.3. NH2 變換為 82A ...................................................................................... 39. 7.

(10) 3.4 Chapter 4. 活化位置環境與氫鍵分析 ....................................................................... 41 結論 (Conclusion)............................................................................... 45. 第二部分分子嵌合模擬 (Part2: Molecular Docking Simulation) ................................................................... 48 Chapter 1. 緒論 (Introduction) ............................................................................ 50. 1.1. 前言 .......................................................................................................... 50. 1.2. 纖維母細胞生長因子受體一與癌症的關聯(The relationship between Fibroblast Growth Factor Receptor 1 and Cancer) ............................... 50. 1.3 Chapter 2 2.1. 研究目標 (Research Objective) .............................................................. 53 模擬方法 (Simulation Methods) ........................................................ 54 分子嵌合模擬 (Molecular Docking)....................................................... 54. 2.1.1. GOLD 與基因演算法 ............................................................................. 54. 2.1.2. 評分函數 ................................................................................................ 56. 2.1.3. 豐富指數與 DUD 資料庫 ...................................................................... 58. 2.2 Chapter 3. 實驗設定 (Experimental Setting) ........................................................... 59 結果與討論 (Results and Discussion) ................................................ 60. 3.1. 側鍊擺動 .................................................................................................. 60. 3.2. 氫鍵位置限制........................................................................................... 63. 3.3. 晶體結構評估........................................................................................... 68. 3.4. 條件收斂性 .............................................................................................. 72. Chapter 4. 結論 (Conclusion)............................................................................... 76. 第三部分總結 (Part3: Concluding Remark) .................................................................................... 78 第四部份附錄與參考資料 (Part4: Appendix and Reference) ............................................................................. 82 Chapter 1. 附錄 (Appendix) ................................................................................. 84. 1.1. 名詞縮寫及翻譯 (Abbreviation and Translation of Nouns) .................. 84. 1.2. 複合物相 TI-2, 所對應之 dv/d及自由能差(G) ................................. 85. 1.2.1. CF3 變換為 82A..................................................................................... 85. 1.2.2. OHW 變換為 82A .................................................................................. 86 8.

(11) 1.2.3 1.3. NH2 變換為 82A .................................................................................... 87 豐富指數分布(Enrichment Factor Distribution) ................................... 88. 1.3.1. 3JS2 A Chain .......................................................................................... 88. 1.3.2. 1AGW .................................................................................................... 89. 1.3.3. 1AGW-Complete Sequence .................................................................... 90. Chapter 2. 參考文獻 (Reference) ......................................................................... 91. 9.

(12) 10.

(13) 圖目錄 (LIST OF FIGURES) Figure 1: Ras/Raf/MEK/ERK 訊息傳導路徑 6 ............................................................ 19 Figure 2: ERK2 蛋白質結構 ....................................................................................... 20 Figure 3: 2OJJ11 蛋白質序列(綠: Beta Sheet；黃: Alpha Helix) ................................ 20 Figure 4: 方法精確度與 CPU 時間之關係圖 ............................................................ 24 Figure 5: 熱力學循環................................................................................................. 25 Figure 6: 溶液(左)及蛋白質中(右)的分子................................................................. 27 Figure 7: 溶液相中 CF3 變換為 82A 熱力學積分..................................................... 33 Figure 8: 複合物相中 CF3 變換為 82A 熱力學積分(TI-1、TI-3)............................. 34 Figure 9: 溶液相中 OHW 變換為 82A 熱力學積分 .................................................. 36 Figure 10: 複合物相中 OHW 變換為 82A 熱力學積分(TI-1、TI-3) ........................ 37 Figure 11: 溶液相中 NH2 變換為 82A 熱力學積分 .................................................. 39 Figure 12: 複合物相中 NH2 變換為 82A 熱力學積分(TI-1、TI-3) .......................... 40 Figure 13: (A)82A、(B)CF3、(C)OHW、(D)NH2 與重要胺基酸之間的三度空間關係 .................................................................................................................. 42 Figure 14: OHW 在 ERK 活性位置中的環境............................................................. 44 Figure 15: 2011 年 Kim, M. H. et. al 等人研究中 26 所建議取代的官能基 ................ 46 Figure 16: FGFR1 在生物細胞訊息傳導路徑所扮演的角色 35 ................................. 51 Figure 17: FGFR1 酪胺酸蛋白質激酶區的 T 結構(pdb code:3JS239) ........................ 52 Figure 18: 3JS2 蛋白質序列(綠: Beta Sheet；黃: Alpha Helix) ................................. 52 Figure 19: 基因演算法之流程圖 45 ............................................................................ 55 Figure 20: Gmetal、Glipo 中 f(r)與 R1、R2 之關係 47 .......................................... 57 Figure 21: FGFR1 蛋白質序列比較:1AGW1 及 3JS2 ................................................. 59 Figure 22: 擺動 3JS2 b-chain 中 10Å 內側鍊擺動與 active ligand 的挑選數量分佈圖 .................................................................................................................. 60 Figure 23: 3JS2 活化中心 Lys514、Ala564 及原結晶結構分子 VM1 的相對位置 ... 62 Figure 24:氫鍵限制後的 active ligands 挑選數量分佈圖 ........................................... 64 Figure 25:擺動 Lys514 及氫鍵限制後的 active ligands 挑選數量分佈圖 .................. 66 Figure 26: GOLD 軟體所再現的小分子結構 ............................................................. 69 11.

(14) Figure 27: SU2 分子在 1AGW 蛋白質結構中的位置與本身之結構缺陷 ................. 69 Figure 28: Top 1%之收斂性分佈圖 ............................................................................ 74 Figure 29: Top 5%之收斂性分佈圖 ............................................................................ 75. 12.

(15) 表目錄 (LIST OF TABLES) Table 1: 完整熱力學循環(Thermodynamics Cycle)模擬的軌跡數量 ........................ 29 Table 2: 82A、33A 分子之間的結構差異 .................................................................. 30 Table 3: 先前研究中所得到的自由能-自由能差....................................................... 30 Table 4: 本研究中類似物的分子結構 ....................................................................... 31 Table 5: 溶液相中 CF3 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G)............. 33 Table 6: 複合物相中 CF3 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G)......... 34 Table 7: CF3 變換為 82A 之自由能-自由能差(G)................................................. 35 Table 8: 溶液相中 OHW 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) .......... 36 Table 9: 複合物相中 OHW 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) ...... 37 Table 10: OHW 變換為 82A 之自由能-自由能差(G) ............................................ 38 Table 11: 溶液相中 NH2 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) .......... 39 Table 12: 複合物相中 NH2 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) ...... 40 Table 13: NH2 變換為 82A 之自由能-自由能差(G) .............................................. 40 Table 14: MD 模擬後所有分子與蛋白質中重要胺基酸間的氫鍵............................. 43 Table 15: 指導教授孫英傑教授與陳博晉同學所變換的分子結構類似物 ............... 45 Table 16: 分子結構變換時溶液相、複合物的自由能差及自由能-自由能差 .......... 46 Table 17: 擺動 3JS2 b-chain 中 10Å 內側鍊擺動與 active ligand 的挑選效率.......... 61 Table 18: 氫鍵限制後的 active ligands 挑選效率 ...................................................... 63 Table 19: 結構擺動 Lys514 及氫鍵限制後的 active ligands 挑選效率 ..................... 65 Table 20: 四個結晶結構設定不同條件的 active ligands 挑選效率 ........................... 67 Table 21: 3JS2 與 1AGW 結晶結構蛋白質的體積 .................................................... 68 Table 22: 各種設定條件對四個結構挑選活性分子的評估(No. of Ligands) ............. 71 Table 23: 各種設定條件對四個結構挑選活性分子的評估(EF) ............................... 71 Table 24: 各種設定條件對四個結構挑選活性分子的評估(EF Ratio to Rigid) ......... 71 Table 25: 3JS2 b-chain 在 Operation=10%時 active ligands 的挑選效率 .................... 73 Table 26: 3JS2 b-chain 在 Operation=10K 時 active ligands 的挑選效率 ................... 73 Table 27: 3JS2 b-chain 在 Operation=30K 時 active ligands 的挑選效率 ................... 73. 13.

(16) 符號(SYMBOLS) V: total potential energy kl: vibrational constant for stretching kθ: vibrational constant for bending Vx: potential energy constant of torsional energy r: distance between atoms q: point charge ε0: dielectric constant ε: a depth of Van Der Waals well σ: Van Der Waals radius F: force m: particle mass ̈ : the second derivative of r under specific time interval , i.e.. ⁄. : Boltzmann constant coupling parameter for TI 1⁄. . H: Hamiltonian energy function in the present system H0: Hamiltonian energy function for initial state H1: Hamiltonian energy function for final state G0: Gibbs energy in ligand 0 binding G1: Gibbs energy in ligand 1 binding Gsoln: Gibbs energy between ligand 0 and ligand 1 in solution phase Gprot: Gibbs energy in Ligand 1 and ligand 1 in complex phase p: linear momentum qposition , ,. Van Der Waals Softcore Potential for =0 Van Der Waals Softcore Potential for =1. 14.

(17) 15.

(18) 第一部分分子動力學模擬 (Part1: Molecular Dynamics Modeling). 16.

(19) 17.

(20) Chapter 1 1.1. 緒論 (Introduction). 前言 (Preface) 在生物化學中，蛋白質激酶(Protein kinases)為轉移酶(Transferase)中的一種，其. 功能是將生物體中具有高能量的 ATP 分子上的磷酸根轉移到移到受體蛋白質上，這種過程稱作磷酸化轉移 (Phosphorylation) ，磷酸化轉移是細胞訊息的傳遞 (Cell signaling)的一部分。. 細胞訊息的失調會導至許多重要的疾病，特別是多種類型的癌症. 3. 以及發炎反. 應。因此，蛋白質激酶在目前為藥物學家重視的癌症藥物標的(Drug target)之一 4。在人類基因體(Genome)中有超過五百個不同的蛋白質激酶存在，這些蛋白質激酶可以將磷酸根接在的幾個不同胺基酸上，分為酪胺酸(Tyrosin)、絲氨酸(Serine)、蘇胺酸(Threonine),分別稱為酪胺酸(Tyrosin)、絲氨酸(Serine)、蘇胺酸(Threonine)蛋白質激酶。雖然這些蛋白質的胺基酸次序不同，但蛋白質激酶的 ATP 結合位置以及整體的三度空間結構卻有著高度的相似性。由於其醫學上的重要性，設計有效抑制蛋白質激酶的小分子藥物便成為藥廠及生技公司當前的目標之一。在這本論文中的第一部分，我們採用電腦輔助藥物設計的方式，推估藥物潛力分子與蛋白質之間的結合強度，並且更進一步設計結合效果優於現有抑制劑的小分子藥物。. 1.2. 細胞外調節激酶二與癌症的關聯 (The Relationship Between Extracellular Signal Regulated Kinase Protein Kinase 2 and Cancer) 在細胞中，Ras/Raf/MEK/ERK訊息傳導路徑接收細胞外的部分訊息，並且將訊. 息傳自細胞核內，此訊息傳導路徑與細胞核中基因表現(Gene expression)、細胞凋零 (Cell apoptosis)、細胞週期(Cell cycles)、分化(Differentiation)、增生(Proliferation)等功能的有關。之中包含了蛋白質激酶(Protein kinase)、轉錄因子(Transcription factor)、細胞凋零調節者(Apoptosis regulator)、凋亡蛋白酶(Caspase execuitioner )等不同功能的蛋白質。這些蛋白質會藉由磷酸化激活(Activate d by phosphorylation)或者失去原 18.

(21) 有的活性(Inactivated by phosphorylation)。5. Figure 1: Ras/Raf/MEK/ERK 訊息傳導路徑 6. 有絲分裂活化蛋白質激酶可分為兩種亞型(Isoform)，分別為有絲分裂活化蛋白質激酶一、有絲分裂活化蛋白質激酶二 (ERK1、ERK2)，這兩種在序列上幾乎有 90% 的相似度。有絲分裂活化蛋白質激酶一 (Mitogen-activated protein kinase 1, MAPK1) 又可被稱作細胞外調節激酶二 (Extracellular signal-regulated kinase 2, ERK2)。在細胞中，ERK 為 Ras/Raf/MEK/ERK 訊息傳遞路徑下游的蛋白質，其功能與基因表現的調節及細胞凋零有著直接的相關。因此在研究中也發現此蛋白質與癌症 7、阿茲海默症 8、侏儒症 9 的生長調節有關。. 19.

(22) ERK2 為一類型的酪胺酸及蘇胺酸蛋白質激酶(Tyrosine and Threonine Kinase )在三度空間 X 光結構上，如 Figure 2 可分為幾個部分: 負責便是結合受體特異性的 ED (residues 157-158) 、 CD (residues 316-319) ；異位結合位 (Allosteric site) 中的 DFG( residues 165-167)、DEF(residues 231-235)；催化片段(Catalytic loop, residues 147-152)；活化片段(Activation loop, residues 171-185)；多甘胺酸片段(Glycine-rich loop, residues 32-36)；絞鍊片段(Hinge, residues 105-109)；MKI 螺旋片段(MKI helices, residues 244-267) 10。這個蛋白質中的 ATP 結合位置約在 Glycine-rich Loop、Hinge、 Catalytic Loop 之間。. Figure 2: ERK2 蛋白質結構. Figure 3: 2OJJ11 蛋白質序列(綠: Beta Sheet；黃: Alpha Helix) 20.

(23) 1.3. 研究目標 (Research Objective) 在先前的研究中，我們已經利用熱力學積分的方法驗證利用此方法對兩個已知. ERK2 抑制劑的 Ki 值評估，證明使用此方法所計算出的自由能-自由能差值(G)優於 MM-PBSA 方法。本篇論文中，我們從原本已知抑制 ERK 較好的小分子開始，進一步設計不同的化學官能基。並且使用分子動力學模擬中熱力學積分的計算精確地計算不同官能基間相對自由能差，探討一系列相似結構之小分子與 ERK 之間結合 (Binding)的差異，並且藉由此方法，設計出比原本抑制劑有更強結合效果之分子。. 21.

(24) 模擬方法 (Simulation Methods). Chapter 2 2.1. 理論. 2.1.1. 分子力學與分子動力學. 現今許多實際化學的問題上，往往系統中的原子數量是非常龐大的。在大尺度的分子系統中，若要使用量子力學方法(Quantum Mechanics)來描述整個系統電子結構(Electronic Structure)則顯得困難且耗時。化學家為了描述此類複雜度高、原子數量極大的系統而發展出分子力學方法(Molecular Mechanics)。分子力學中，我們不考慮電子的運動行為，我們僅考慮以原子核為主所呈現出的特性。在分子力學計算裡面，我們使用古典的力場模型 (Force Field) 取代量子力學中的薛丁格方程式 (Schrödinger equation)來描述系統的能量。Eq.1 為典型古典力場描述系統中分子位能的函數。. V( ⃑) = ∑ + ∑. ( − 1+. ,. (. ) + ∑ − ) + ∑. (. − ∑. ,. (. ) +4. [(. ). −(. ) ]). Eq. 1 在此函數中，我們將系統的位能分為：1.以虎克定律描述兩原子簡諧振子的振動位能(Stretching) 2.同樣以簡諧振子方式描述三個原子間的彎曲位能(Bending) 3.描述四個原子間兩面角的位能(Torsional Energy) 4.以庫倫位能(Coulomb Potential)及蘭納 - 瓊斯位能 (Lennard-Jones Potential) 分別描述距離三個化學鍵以上的靜電 (Electrostatic)及凡得瓦作用力(Van Der Waals’ interaction)。. 有了力場函數之後，我們便可以得知整個系統的位能高低。由於分子結晶結構的原子座標並非完全符合力場模型中位能曲面(Potential Energy Surface)能量低點的結構，因此能量極小化(Energy Minimization)則被引進計算流程中。藉由此程序我們可以得到結構中力場位能的區域低點(Local Minima)，使得分子結構相對符合力場模型。在此我們所採用的是最陡下降法(Steep Decent Method)12 與共軛梯度法(Conjugate Gradient Method) 13 來實行此程序。. 22.

(25) 藉由電腦計算在不同時間點時，系統中原子的運動行為，這種方法稱作分子動力學(Molecular Dynamics)。分子動力學模擬中，我們將得到隨著時間變化時，用來描述粒子位置及速度的軌跡。我們的系統裡考慮的是整個蛋白質大分子，此類大分子中量子效應相對顯得微小，故我們可以使用牛頓第二定律來描述粒子的行為。. 牛頓方程式(Newton’s Law)中描述的是隨著不同作用力所產生連續的運動行為，其函數如 Eq.2 : (⃑ ) =. ̈ =−. ( ⃑). Eq. 2 有了力場描述系統的位能後，再將此位能對位置微分即可得知瞬間的系統作用力。在一剛開始解牛頓方程式的過程中，我們必須利用亂數產生波茲曼分佈 (Boltzmann Distribution)以得到特定溫度下的初始速度(Initial Velocity)。有了位能函數、初始位置及初始速度，我們將得知在不同時間點下系統的變化，並且得到一連串用來描述系統整體行為的軌跡(Trajectories)。從這些軌跡中，便可得知在 Ergodic Hypothesis 底下真實平衡系統中的溫度、壓力、密度、相對位能等物理量。. 23.

(26) 2.1.2. 自由能計算. 要如何準確計算分子間結合的自由能(Gibbs Free Energy)及溶劑效應一直是複雜系統的計算化學領域所面臨的挑戰。在過去的二十年間計算化學家也發展了許多方法來做這一方面的計算。在蛋白質與小分子結合(Binding)的問題中，這些方法包括 Docking14、3D-QSAR15、MM-GB/PB SA16、FEP17 及 TI18,19，這些方法根據原理的不同而估計出來的自由能大小也會不同。一般而言，越接近統計熱力學上第一原理(First Principle in Statistical Thermodynamics)以及假設越少的方法所得到的自由能越接近真實情況；相對的，這些越準確的方法所耗費的 CPU 時間也更為可觀。Figure 4 為各種方法精確度與計算 CPU 時間之示意圖。. Figure 4: 方法精確度與 CPU 時間之關係圖. 24.

(27) 2.1.3. 熱力學循環與熱力學積分. 在熱力學(Thermodynamics)中，自由能為狀態函數(State Function)。因此我們欲得知物體兩狀態之間的自由能差，則可藉由計算不同的路徑而得到兩狀態間相同的自由能差值。把這個概念引導到蛋白質與小分子之間的相互結合情形，如 Figure 5 所示:. Figure 5: 熱力學循環. 在 Figure 5 中，生物實驗可測得知數值為不同抑制劑與蛋白質結合時的 Ki 值，從 Ki 值中我們可將其透過 Eq. 3 轉換為自由能差(G0、G1). G = -RTln(1/Ki) Eq. 3. 然而在計算化學中，若要直接計算分子與蛋白質之間的自由能差(G0、G1)顯得極度困難，因此化學家計算Gsoln、Gprot。藉由計算Gsoln、Gprot 的差值G 亦可得知G0、G1 間的差值。(G=G1-G0=Gprot-Gsoln)然而真實物理情況下 Gsoln、Gprot 並不會發生，所以這部份的計算又被稱作 Alchemical Calculation，計算的方法包含 TI、FEP。. 在統計力學中(Statistical Mechanics)，配分函數(Partition Function) 可用來所描述的為封閉系統中(Canonical Ensemble)特定溫度時的各種能階狀態的機率分布情況， 25.

(28) 所有物理量如總能(Total Energy)、自由能、焓(Enthalpy)、熵(Entropy)皆可從配分函數(Partition Function)得知。古典的配分函數可由 Eq. 4 表示 Q = Eq. 4 而平衡狀態的亥姆霍茲自由能(Helmholtz free energy)則可以表示為: =−. １. lnQ. Eq. 5 而熱力學積分便是奠基在此表示式之下。在熱力學積分法. 18. 中，我們將描述系. 統總能的 Hamiltonian Function, H 表示為下列函數:. H() = (1-)H0+H1 Eq. 6. 從 Eq.6 中我們可以發現，當為 0 時，我們所描述的是初始狀態的 Hamiltonian Function；反之當為 1 時，我們所描述的是最終狀態的 Hamiltonian Function。描述完位能形式後，引入公式 Eq.4、Eq.5，我們可知自由能的變化(A)可表示為下式 Eq.7:  ∆A =. ( ). ( )=. ⟨. =. ( ). ⟩. 及 ∆G =. ⟨. ( ). ⟩. =. ⟨. ( ). ⟩. Eq. 7. 因此我們可以將等於 0 到 1 之間的所有平衡狀態值所對應到的. ( )⁄. 加總，. 則可得知從=0 到=1 之間變換過程的G soln、Gprot 以及G 值。得知G 後，我們即可從 Eq. 3 得知轉換後分子與原本具有實驗值之參考分子之相對 Ki 值，便可以應用此計算方法推估未知化合物與蛋白質結合時的平衡常數大小，以達到節省合成小分子及酵素實驗所需之金錢與時間。. 26.

(29) 模型設定 (Modeling Setup). 2.2. 由於本篇在熱力學積分類型的分子動力學模擬計算中，我們分別需要建置 1.在溶液中的分子 2.在蛋白質複合體中的分子。以下分為兩個小節敘述。. 2.2.1. 溶液的分子. 從蛋白質資料庫(RCSB，Protein Databank)下載 82A 分子在蛋白質中 X 光結晶結構(pdb code: 2OJJ)中的 3D 化學結構，利用 Ds Viewer Pro 將原本 82A 的官能基改變成其他欲變化的化學結構，再利用 Gaussian03/HF/6-31G*計算此類小分子抑制劑的電荷分佈。接著透過 AMBER11 Tool 中的 tLeap 模組將此類分子模型浸至 41.5Å x 28.7 Å x 39.3 Å 的邊界水盒子(Periodical Boundary Condition Solvent box)中，其中包含了 1451 個 TIP3P 水分子，整個系統共包含約 4401 個原子。並且產生 AMBER11 程式所需要的計算格式。. 2.2.2. 蛋白質中的分子. 有了各式小分子的電荷分佈後，我們將各種不同的小分子放入蛋白質原始 X 光結晶結構(pdb code: 2OJJ)中。並且在 tLeap 中，我們將蛋白質結構內所欠缺的氫原子依照標準的鍵角鍵長補足。接著將整個系統浸入 92.2 Å x 74.3 Å x 97.7 Å 的水盒子中，並且加入適當的 Na+離子使得整個溶液系統呈現電中性。這其中包含了 16322 個 TIP3P 水分子，整個系統共包含約 54661 個原子。. Figure 6: 溶液(左)及蛋白質中(右)的分子 27.

(30) 2.2.3. 模擬參數. 此篇分子動力學計算中，我們所使用的軟體為 AMBER11/Sander20 程式集。在描述力場上，蛋白質部份我們所使用的是 AMBER ff03.r121 力場；小分子所使用的是 GAFF(General AMBER Force Field) 力場來描述粒子間的作用力。. 由於標準的熱力學積分程序為了克服收斂性及端點歧異性(Endpoint Singularity) 所以必須將轉換程序分為三個部分 :TI-1(Charge off) 、 TI-2(Van Der Waals Transformation)、TI-3(Charge on)。每一個部分中，我們計算、、、 共九條軌跡，並且從每一條軌跡中得到一個 dv/d值。值得一提的是，在 TI-2 的部分，由於牽涉原子的消失與出現，所以整體位能的描述也必須特別小心。在這邊，我們為了處理這一類問題，在這一部分的凡得瓦位能描述上，我們所使用的為 Soft-Core Van Der Waals Potential22。其位能函數如 Eq. 8:. = 4 (1 − ). ,. ⎧ ⎪ ⎨ ⎪ ⎩. 1. −. 1 +. +. ⎫ ⎪  ⎬ ⎪ ⎭. and. ,. =4. ⎧ ⎪ ⎨ ⎪ ⎩. 1 (1 − ) +. −. 1 (1 − ) +. ⎫ ⎪  ⎬ ⎪ ⎭. Eq. 8  TI-1、TI-3 的部分我們使用 SHAKE23 演算法，因此時間梯步為 2fs；而 TI-2 因為要做凡得瓦轉換而無法使用 SHAKE 演算法固定氫原子的鍵長，因此將時間梯步改為 1fs。. 無論在 TI-1、TI-2、TI-3 中，我們先使用最陡下降法、共軛梯度法做能量極小化，再進行 MD 模擬。在 MD 模擬中 1.使用 Langevin thermostat24 將系統做 10K 升溫到 300K 的 NVT ensemble 模擬 10ps，以防止氣泡(Bubble)產生，再將系統維持在 300K 持續 50ps。2.做完 NVT ensemble 之後，我們改做 300K、1atm 的 NPT ensemble 模擬 200ps 作為 equilibration run，最後做兩段分別為 600ps 的 300K、1atm 的 NPT 28.

(31) ensemble 作為收集數據點的 production run (但水中的分子由於環境較蛋白質中的分子單純，只做一次 600ps 的 production run)。在整個過程中，我們使用 PME 估計長程的靜電作用力(Electrostatic interaction)，並且將 Van Der Waals Cutoff Distance 設為 9 Å。無論溶液相(Solution Phase)或者複合物相(Complex Phase)皆是取最後 600ps 的 production run 上的樣本數據點(Sampling)作為 TI 的平均。. 蛋白質中小分子周圍的環境相對水中來的複雜，所以 TI-2 凡得瓦轉換計算比較不容易得到收斂，因此 TI-2 的部分我們分別用六個不同的亂數種子(Random number seed)跑六條軌跡，並且將此六條軌跡所分別積分後之面積平均值作為此部分的G。  Table 1: 完整熱力學循環(Thermodynamics Cycle)模擬的軌跡數量 Condition. TI-1. TI-2. TI-3. Ligand in Solution. 9x2. 9x2. 9x2. Ligand bounded Protein. 9x2. 9x6x2. 9x2. 由 Table 1 可以知道在我們的熱力學積分模擬研究中，所使用的 CPU 時間是非常的龐大，我們所使用的電腦叢集為台灣高速網路與計算中心的 Irish-IBM1350 clusters。註：在此所提到的 X 分子分別為 Chapter3 所提到的 CF3、OHW、NH2 三種類型的分子結構。. 29.

(32) Chapter 3. 結果與討論 (Results and Discussions). 2007 年 Alex M. Aronov 在文獻 11 中利用生化實驗的方式測出了一系列 ERK 蛋白質的抑制劑的 Ki 值並在同篇文獻中作者亦提供了兩個結晶結構(2OJJ、2OJI)。另外 2009 年 Alberto Del Rio 等人利用 MM-PBSA 的方式 25，計算出一系列分子的自由能。在我們實驗室先前的研究中，我們利用熱力學積分的方式，針對兩個不同蛋白質-抑制劑(Protein-Ligand Complex)結晶結構(82A in 2OJJ、33A in 2OJI)，計算此兩個抑制劑之間的自由能-自由能差(G)。. Table 2: 82A、33A 分子之間的結構差異 H N. O. H N. HN. N. Cl. 82A. 33A. 發現我們利用熱力學積分的方式所計算出來的自由能-自由能差比起 Alberto Del Rio 利用 PBSA 方式算出來的自由能-自由能差更能接近 2007 年 Alex M. Aronov 測出 Ki 值轉換後的自由能-自由能差。如 Table 3 所示:. Table 3: 先前研究中所得到的自由能-自由能差引用文獻. 自由能-自由能差(G), Kcal/mol. 2007, Alex M. Aronov. -2.24 (從 Ki 實驗所轉換的G). 2009, Alberto Del Rio. -5.65. 我們先前的研究. -1.6. 30.

(33) 從 Table 3 中我們可以知道，在 2009 年 Alberto Del Rio25 所做的研究中，模擬與實驗轉換後的自由能-自由能差的差值約為 3.4 Kcal/mol；在我們先前的研究中，我們所算出來的自由能-自由能差為-1.6 Kcal/mol，與 2007 年生化實驗中從 Ki 值所轉換的自由能-自由能差的差值約為 0.6 Kcal/mol。在藥物設計、蛋白質與小分子結合 (Protein-Ligand Binding)的領域中，G 的差值若為 1 Kcal/mol，則對應到原本 Ki 值約為十倍(Eq.3)。因此準確的評估自由能的大小對藥物設計而言為一個亟需改善的課題。從上面的數據中，不難看到我們所使用的熱力學積分方法在評估自由能上優於 MM-PBSA 方法。在此篇論文中，我們站在實驗室先前的研究上，進一步利用熱力學積分的方式，精確地設計結合效果優於 2007 年在 Alex M. Aronov 所發表的文獻中未提到的分子架構(Molecular Structure)，並且以此方法為開端，希望能夠利用熱力學積分方法準確的預測化合物，以減少有機化學家在探索藥物及合成時所需要大量的金錢與時間。在此列出所有我們轉換的分子結構:. Table 4: 本研究中類似物的分子結構. Cmpd. R. Cmpd. 82A. CF3. OHW. NH2. 31. R.

(34) 在正常的熱力學積分模擬中，起始分子( =0)設定為原始已知結晶結構的分子，而終點分子( =1)設定想預測的分子。但是在我們的研究中，為了防止 TI 產生大取代基時，蛋白質結構中的原子與小分子重疊而造成位能太高的情形，我們將比較大的分子 CF3 作為起始分子( =0)，而將結構較小的 82A 設定為終點分子( =1)。在以下三個部分: CF3 轉換為 82A、OHW 轉換為 82A、NH2 轉換為 82A 皆以這種方式進行。由於技術上的細節，導致轉換方式與直觀物理化學的變化相反，所以在熱力學循環中所計算出來的G 會與直觀的方式相反。也就是說，本篇論文中若G 值大於零，則代表新設計出的化合物與蛋白質間結合的情形優於原本結晶結構中 82A 分子；反之，若G 值小於零則代表新設計出化合物之結合能(Binding Energy)劣於 82A。. 32.

(35) CF3 變換為 82A. 3.1. 第一個部分我們將我們所設計的官能基 CF3 透過熱力學積分程序變換為原本結晶結構 2OJJ 中的 82A 分子，在文獻中 1182A 為目前對 ERK2 最好的抑制劑，其平衡常數 Ki 值為 2nM。其在溶液相中的自由能差如 Table 5。從 Table 5 中的 Total 部份我們可以明顯的發現，單獨在水溶液條件下 82A 比起 CF3 更容易存在水中，其中的原因可能是 82A 中的 F、Cl 官能基比起 CF3 官能基更親水，這其實很符合我們在物理化學上面的直覺。其水溶液相中 dv/d對的原圖如 Figure 7 所示。從 Figure 7 中可以發現 TI-1、TI-3 兩條單純移除電荷的條件比 TI-2 做凡德瓦轉換的曲線平緩，從這個地方可以得知做凡德瓦轉換時，dv/d與值並沒有特別的相關性。. Table 5: 溶液相中 CF3 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) Step. TI-1. TI-2. TI-3. Total. Gibbs Free Energy (Kcal/mol). -18.42. 17.98. -5.34. -5.78. CF3 to 82A in solution phase 30 20. dv/d. 10 0. TI-1. -10. TI-2 TI-3. -20 -30 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9.  Figure 7: 溶液相中 CF3 變換為 82A 熱力學積分. 33. 1.

(36) 在這裡所要描述的是 CF3 變換為 82A 在複合物環境下的變換情形。如同水溶液的狀態，Table 6 中 CF3 變換為 82A 的自由能差依然為負值，同樣代表 82A 較 CF3 更適合存在 ERK 活化位置中，也就是說 ERK 的結晶位置可能有著與 82A 相同的環境(例如: 親疏水性、空間障礙度、電荷分布情形等等)。並且從 Figure 8 中發現， TI-1、TI-3 的曲線不像水溶液中的 TI-1、TI-3 平緩，在這裡可能的原因是小分子在蛋白質活化位置中的環境遠比水溶液中的環境複雜許多，因此在不同軌跡中的變化幅度也有明顯的特異性存在。. Table 6: 複合物相中 CF3 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) Step. TI-1. TI-2. TI-3. Total. Gibbs Free Energy (Kcal/mol). -20.33. 23.29. -9.27. -6.31. (註:關於此 TI-2 為六條軌跡的平均資料，詳見本論文第四部份附錄 1.2.1). CF3 to 82A in complex phase 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 1. dv/d. 0 -10 TI1. -20. TI3 -30 -40.  Figure 8: 複合物相中 CF3 變換為 82A 熱力學積分(TI-1、TI-3). 分別有了水溶液中及複合物中的G 之後，我們進一步可以得知的是兩個狀態之間的自由能-自由能差 G，根據 2.1 節實驗理論，我們從G 可以得知在做熱力學積分研究時，變換前與變換後哪一個分子與蛋白質的結合能力較佳，其結果如 Table 7:. 34.

(37) Table 7: CF3 變換為 82A 之自由能-自由能差(G). Gibbs Free Energy (Kcal/mol). Solution G. Complex G. G. -5.78. -6.31. -0.53. 從 Table 7 中我們可以知道，這邊的G 為負值，所以代表 82A 比 CF3 更適合此蛋白質環境，至於什麼原因造成這種現象，將會在 3.4 節做更進一步的結構分析。. 35.

(38) OHW 變換為 82A. 3.2. 在 OHW 變換為原本結晶結構中分子 82A 的這水溶液部分，從 Table 8 中看的出來 TI-1 的G 為正值，這代表我們將 OHW 移除電荷時對整個系統來說是不趨向這種狀況的。從化學上來說，因為 OHW 上面的氫氧基與整體水溶液環境較為接近，兩者皆為具有極性的化學結構；在 TI-2 的部分發現凡德瓦轉換時，中性 OH 官能基與水之間的作用力型態亦比 82A 中中性的 Cl、F 原子更適合存在於水中；而 TI-3 環境下我們可以發現，加回電荷後 82A 分子的 Cl、F 原子比起不帶電的 Cl、F 原子更適合水環境。. 另外從 dv/d對圖(Figure 9)中也看得出來，如同 3.1 節中的 Figure 7，我們可以觀察到 TI-2 曲線的不如 TI-1、TI-3 來的平緩，再次驗證了凡德瓦轉換過程中，每一條軌跡的變動性比起電荷移除(TI-1)或添加電荷(TI-3)來的大。. Table 8: 溶液相中 OHW 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) Step. TI-1. TI-2. TI-3. Total. Gibbs Free Energy (Kcal/mol). 26.24. 3.57. -5.34. 24.47. OHW to 82A in solution phase 40. dv/d. 30 20 10. TI-1. 0. TI-2. -10. TI-3. -20 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9.  Figure 9: 溶液相中 OHW 變換為 82A 熱力學積分 36. 1.

(39) Table 9 中所說的是在 ERK2 蛋白質結晶活化位置中 OHW 變換為 82A 時的自由能差。比較 Table 8、Table 9 之後我們可以看的出來在 TI-1 條件下，蛋白質活化位置的環境比起水環境更不趨向讓原本有帶電的 OH 官能基變為沒有帶電的 OH 官能基，從數據上說明了活化位置的環境比起水環境來的更親水；而在 TI-2 的部分，比起 Table 8 的 TI-2 我們可以發現在單純凡德瓦轉換，由於蛋白質中的變換牽涉到環境中的立體障礙，所以複合物中入要將 OHW 變換為 82A 原子更顯得困難；比較 Table 8、Table 9 的 TI-3 部分，由於此蛋白質活化位置的環境為親水環境並且附近的官能基比起水來更為固定、緊密，因此造成複合物環境下 TI-3 的自由能差更為負值。. 從 Figure 10 及本論文第四部份 1.2.2 節中的 dv/d圖看來，我們可以明顯得知無論在 TI-1、TI-2 及 TI-3 中的曲線，都不像溶液中的曲線來的平順，其中最大的原因是蛋白質環境遠比水溶液環境來的複雜，每條軌跡間的座標相異程度也更大，所以造成積分曲線的顛簸。. Table 9: 複合物相中 OHW 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) Step. TI-1. TI-2. TI-3. Total. Gibbs Free Energy (Kcal/mol). 28.21. 7.61. -9.27. 26.55. (註:關於此 TI-2 為六條軌跡的平均資料，詳見本論文第四部份附錄 1.3.2). OHW to 82A in complex phase 40. dv/d. 30 20 10 TI1. 0. TI3. -10 -20 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 1.  Figure 10: 複合物相中 OHW 變換為 82A 熱力學積分(TI-1、TI-3) 37.

(40) 最後我們將複合物相的自由能差減去溶液相的自由能差，如 Table 10。與 3.1 節 CF3 不同的是，這一節中我們所得到的G 為正值。這代表在 ERK 中，我們新設計的 OHW 分子比原本結晶結構中的 82A 分子與蛋白質有著更強的結合能力。我們也發現若將苯環上面的官能基改為極性類官能基 OH 其顯示的自由能與替換成 CF3 官能基有著不同的效應。. Table 10: OHW 變換為 82A 之自由能-自由能差(G). Gibbs Free Energy (Kcal/mol). Solution G. Complex G. G. 24.47. 26.55. 2.08. 38.

(41) NH2 變換為 82A. 3.3. 由於 3.2 節的結果告訴我們，若將苯環上面的官能基取代為極性官能基，則可以增加蛋白質 ERK 與小分子的結合強度。所以在 3.3 節中，我們所變換的分子是從 NH2 到 82A，把原本在 3.2 節所設計的 OH 官能基改為 HN2 官能基。從 Table 11 中可以看出當 NH2 變換為 82A 在 TI-1(移除 NH2 電荷)時，與 OHW 相同的是這一步的自由能差為正值，符合極性官能基變為中性時的比較不趨向中性原子的直覺；TI2、 TI3 兩個部分和 OHW 變為 82A 分子的情形類似，所以就不在這裡做進一步的討論。不單如此，從 Figure 11 中也可以發現與 3.2 節中的 Figure 9 有類似的效應。. Table 11: 溶液相中 NH2 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) Step. TI-1. TI-2. TI-3. Total. Gibbs Free Energy (Kcal/mol). 40.00. 6.60. -5.34. 41.26. NH2 to 82A in solution phase 50. dv/d. 40 30 20. TI-1. 10. TI-2. 0. TI-3. -10 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 1.  Figure 11: 溶液相中 NH2 變換為 82A 熱力學積分. 在複合物相的部分，比較 Table 11 與 Table 12 也可以發現如同 3.2 節中的想法:ERK 蛋白質活化位置環境在親水性以及結構方面，比水溶液環境下來的親水且結構固定，而造成複合物相中的G 較純溶液中之絕對值來的大。 39.

(42) 不僅如此，Figure 12 及中本論文第四部份 1.2.3 節中的 dv/d圖看來，比起水溶液條件，蛋白質條件中的軌跡更因為環境因數使得曲線較為不平。. Table 12: 複合物相中 NH2 變換為 82A 之 TI-1、TI-2、TI-3 之自由能差(G) Step. TI-1. TI-2. TI-3. Total. Gibbs Free Energy (Kcal/mol). 41.58. 10.32. -9.27. 42.63. (註:關於此 TI-2 為六條軌跡的平均資料，詳見本論文第四部份附錄 1.2.3). NH2 to 82A in complex phase 50. dv/d. 40 30 20 10. TI1. 0. TI3. -10 -20 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 1.  Figure 12: 複合物相中 NH2 變換為 82A 熱力學積分(TI-1、TI-3). 有了複合物相及溶液相的G 之後，便可求得G。我們發現當 NH2 變換成 82A 之後，如 Table 13 其G 為正值。這符合我們先前的預測將苯環上面的官能基取代為極性更高的官能基會得到較好的結合能力有著相輔相成的結論，但是其與蛋白質的結合效果並未比 OHW 來的強，關於這點的討論將在下一節中做結構上的說明。. Table 13: NH2 變換為 82A 之自由能-自由能差(G). Gibbs Free Energy (Kcal/mol). Solution G. Complex G. G. 41.26. 42.63. 1.37. 40.

(43) 3.4. 活化位置環境與氫鍵分析在這一節中，我們將探討以上具有不同官能基的四個分子與蛋白質活化位置附. 近胺基酸環境間的關係。四個分子包括:原本 ERK 結晶結構 2OJJ 中的 82A 分子，另外設計 CF3、OHW 及 NH2 分子。在這裡我們將原本用來做 TI 的四個蛋白質複合物檔案分別跑四條單獨 2ns 的軌跡，並且個別分析這四條軌跡中抑制劑分子重要胺基酸殘基(Residue)間的氫鍵比例及相對空間位置。. 從氫鍵的分析 Table 14、Figure13 可以發現，原本結晶結構上 82A 分子已經有很多條顯著的氫鍵。Met106、Gln103、Lys52、Ser151 皆與 82A 有著氫鍵。這些氫鍵也是造成 82A 原本的 Ki 值就已經達到 2nM 的主因。由於這些氫鍵原本已經存在在分子的主體架構中，所以大部分的分子架構不需改變便可得到不錯抑制 ERK2 的效果，因此我們必須進一步設計後續官能基。. 在 CF3 的部分，我們並看不到除了原本 Met106、Gln103、Lys52 三個胺基酸以外的氫鍵，並且原本在 Ser151 上面的氫鍵因為蛋白質結構調整的關係而消失。這可能是 CF3 結合能力較差的緣故。. OHW 的部分發現除了 Met106、Gln103、Lys52 三個胺基酸上面原本的氫鍵之外，還可以發現 OH 官能基會與 Gln69、Tyr34、Lys52 三個官能基產生比例不小的氫鍵作用力，從這三個氫鍵作用力中得知在我們替換的苯環附近其實是以親水性的氨基酸序列為主，由於結構上的氫鍵作用力以其極性高的環境，使得 OHW 與 ERK2 的結合能力比 82A 高約 2 Kcal/mol。雖然化學上面 2 Kcal/mol 是非常小的能量(碳-碳鍵鍵能約為 100 Kcal/mol)，但由於在蛋白質與抑制劑結合的領域中，一千卡的能量就足以讓平衡常數 Ki 質差到十倍，因此對藥物設計而言，結合能力若比起原來已知抑制劑的結合自由能穩定兩千卡，則代表抑制效果為原本抑制劑的一百倍。. 41.

(44) 而 NH2 的結合能力比起 82A 來說更適合作為抑制劑，但比起 OHW 來說結合效果並沒有那麼好。從 Table 14 中發現，與 NH2 官能基形成的氫鍵的胺基酸與 OHW 近乎類似，但是整體來說氫鍵的比例卻沒有 OHW 那麼高。另外一個原因可能是因為 OH 的原子大小、官能基形狀及部分電荷分布所影響的空間效應(Steric effect)及 dipole moment、quadrupole moment 的影響，使得 OHW 的結合效果約比 NH2 來的好 0.8 Kcal/mol 左右。. Figure 13: (A)82A、(B)CF3、(C)OHW、(D)NH2 與重要胺基酸之間的三度空間關係. 42.

(45) Table 14: MD 模擬後所有分子與蛋白質中重要胺基酸間的氫鍵 DONOR res@atom. ACCEPTORH res@atom. ACCEPTOR res@atom. %occupied. distance. angle. 82A@H2 82A@H15 Lys52@HZ3 Lys52@HZ1 Lys52@HZ2 82A@H9. 82A@N1 82A@N3 Lys52@NZ Lys52@NZ Lys52@NZ 82A@O1. 89.05 98.05 30.20 28.55 26.50 66.80. 2.875 ( 0.13) 2.932 ( 0.16) 2.889 ( 0.15) 2.889 ( 0.16) 2.881 ( 0.14) 2.797 ( 0.19). 38.88 (12.98) 17.56 ( 9.15) 29.08 (15.07) 28.30 (14.36) 29.54 (14.86) 18.62 (10.62). CF3@H2 CF3@H15 Lys52@HZ2 Lys52@HZ3 Lys52@HZ1. CF3@N1 CF3@N3 Lys52@NZ Lys52@NZ Lys52@NZ. 95.65 96.45 47.30 21.10 20.70. 2.925 ( 0.16) 2.989 ( 0.18) 2.809 ( 0.11) 2.812 ( 0.11) 2.794 ( 0.11). 22.29 (10.75) 14.51 ( 7.70) 30.35 (13.63) 31.57 (14.01) 32.50 (14.14). OHW@H2 OHW@H15 Lys52@HZ1 Lys52@HZ3 Lys52@HZ2 Lys52@HZ3 Lys52@HZ1 Lys52@HZ2 OHW@H29 Tyr34@HH OHW@H29. OHW@N1 OHW@N3 Lys52@NZ Lys52@NZ Lys52@NZ Lys52@NZ Lys52@NZ Lys52@NZ OHW@O30 Tyr34@OH OHW@O30. 96.65 95.10 33.80 29.80 25.60 10.30 8.50 7.85 89.15 48.15 25.65. 2.870 ( 0.13) 2.964 ( 0.17) 2.791 ( 0.11) 2.776 ( 0.10) 2.785 ( 0.10) 3.311 ( 0.14) 3.307 ( 0.15) 3.292 ( 0.16) 2.638 ( 0.12 ) 3.126 ( 0.21) 3.010 ( 0.36). 27.68 (14.47) 16.41 ( 8.69) 29.36 (14.40) 27.91 (13.85) 28.98 (14.98) 28.67 (15.16) 26.78 (14.00) 27.22 (13.87) 13.81 ( 7.86) 41.42 (11.73) 35.92 (18.76). NH2@H2 NH2@H15 Lys52@HZ1 Lys52@HZ3 Lys52@HZ2 Lys52@HZ2 Lys52@HZ1 Lys52@HZ3 NH2@H49 NH2@H50 NH2@H49. NH2@N1 NH2@N3 Lys52@NZ Lys52@NZ Lys52@NZ Lys52@NZ Lys52@NZ Lys52@NZ NH2@N30 NH2@N30 NH2@N30. 100.00 92.00 57.50 21.50 16.00 8.00 3.00 37.00 55.00 26.00 80.50. 2.877 ( 0.13) 2.967 ( 0.17) 2.859 ( 0.14) 2.832 ( 0.11) 2.790 ( 0.09) 3.295 ( 0.12) 3.176 ( 0.16) 3.122 ( 0.17) 3.061 ( 0.23) 3.111 ( 0.20) 2.903 ( 0.16). 25.60 (12.80) 15.70 ( 8.52) 28.18 (14.07) 26.34 (14.94) 21.80 (10.34) 22.25 ( 9.70) 26.08 (15.69) 28.86 (14.66) 32.93 (15.72) 37.48 (14.98) 24.29 (14.34). 82A Met106@O Gln103@OE1 82A@O2 82A@O2 82A@O2 Ser151@O CF3 Met106@O Gln103@OE1 CF3@O2 CF3@O2 CF3@O2 OHW Met106@O Gln103@OE1 OHW@O2 OHW@O2 OHW@O2 OHW@O30 OHW@O30 OHW@O30 Gln69@OE1 OHW@O30 Gln69@OE2 NH2 Met106@O Gln103@OE1 NH2@O2 NH2@O2 NH2@O2 NH2@N30 NH2@N30 NH2@N30 Gln69@OE2 Tyr34@OH Gln69@OE1. 43.

(46) 此外從 Figure14 的環境圖中我們看得出來:(A)圖中所顯示的是蛋白質表面方式表示 OHW 附近的環境。藍色所代表的為親水性的氨基酸分子，白色為中性或小胺基酸而紅色表示疏水的氨基酸類別。明顯地看的出來在 ERK2 結晶結構活化位置中，大部分胺基酸為親水性胺基酸，特別是我們所取代苯環附近官能基的環境也是以親水性胺基酸為主體。(B)圖為(A)圖以彩帶模式顯示的重要胺基酸位置，更進一步可以發現附近胺基酸來自不同的二級結構(Lys52 in -sheet、Tyr32 in Glycine-rich loop、 Gln69 in -helix)。並且從(A)、(B)兩個圖更能清楚表達活化位置的環境以及重要胺基酸位置的相對位置。. Figure 14: OHW 在 ERK 活性位置中的環境 (A)表面形式表示、(B)彩帶形式表示 44.

(47) Chapter 4. 結論 (Conclusion). 在我們的 TI-MD 研究中，我們以原本 82A 為起始分子，取代了上面的官能基。發現在原本的分子架構下，若將苯環上面間位(meta position)的官能基取代為適當的 OH、NH2 等極性官能基，則可以增加此分子與 ERK2 蛋白質的結合能力(Binding Affinity)。並且，從不同取代基上面氫鍵的分布情形以及蛋白質活化位置的親水、疏水胺基酸殘基的環境分析為何會造成這種結果。除了我的研究部分，根據指導教授孫英傑教授與研究室的陳博晉同學所做的計算指出，若取代後的官能基太大亦不利於其分子與 ERK 蛋白質的結合。Table 15 為 CH3 及 CFO 分子。CH3 與本篇論文中所提到的 CF3 分子有著相似類型皆為疏水性的官能基；CFO 與 OHW 的 OH 取代基位置相同，只是又在對位(para position)新增了立體障礙較大的 CF3 官能基。. Table 15: 指導教授孫英傑教授與陳博晉同學所變換的分子結構類似物. Cmpd. R. Cmpd. CH3. CFO. 45. R.

(48) 根據 Table 16 比較 CH3、CF3:從這兩個皆為非極性類的取代基，雖然兩者的結合效果都不如 82A好，但可以發現比較小的 CH3 基團的G 值稍微好過 CF3 之G；比較 OHW、CFO 兩個分子也可以發現同樣的情況:若單只有 OH 官能基，則G 約為 2 Kcal/mol，如果在對位加上 CF3 官能基則G 只剩下 1 Kcal/mol。從這兩組官能基比較出空間效應對整個官能基的設計也是很重要的因素之一。有趣的是在 2011 年 Kim, M. H. et. al 等人的 3D-QSAR、Docking 研究中 26 指出這一個區域設計的官能基應以較大的官能基為主。他們的研究結果與我們的結果有著根本上的差異，這再次顯示了在針對蛋白質複雜系統時，若用不同的研究方法探討相同問題則可能得到不同的結果。而我們在本篇論文中所使用的熱力學積分方法更能得到接近真實自由能間的差值並且做更仔細的能量分析。. Table 16: 分子結構變換時溶液相、複合物的自由能差及自由能-自由能差 Transformation. Solution. Complex. CF3 to 82A OHW to 82A NH2 to 82A CH3 to 82A CFO to 82A. -5.78 24.47 41.26. -6.31 26.55 42.63. -3.59. -2.52. G -0.53 2.08 1.37 -0.1 1.06. 註:在我們的研究中，若G 大於零越多則代表新設計的分子結構與 ERK2 的結合效果越優於原本結晶結構中的 82 分子；反之亦然。. Figure 15: 2011 年 Kim, M. H. et. al 等人研究中 26 所建議取代的官能基. 46.

(49) 47.

(50) 第二部分分子嵌合模擬 (Part2: Molecular Docking Simulation). 48.

(51) 49.

(52) Chapter 1 1.1. 緒論 (Introduction). 前言在藥物化學領域中，如何有效率並且廣泛地找出各種具有抑制蛋白質、DNA 的. 小分子是一件極困難的事情。在 80 年代末期，為了有效率地找出先導化合物(Lead Compound)，電腦輔助藥物設計中的分子嵌合方法因而被發明。在這一部分的研究中，我們使用這種方法，討論 University of California, San Francisco 所提供的 DUD database 在設定纖維母細胞生長因子受體蛋白質激酶條件時，如何有效率的挑選出具有生物活性的小分子抑制劑。. 1.2. 纖維母細胞生長因子受體一與癌症的關聯 (The relationship between Fibroblast Growth Factor Receptor 1 and Cancer) 纖維母細胞生長因子受體總共包含了四個家族(FGFR 1-4)。FGFRs 在目前的研. 究中發現此類蛋白質與許多生物過程有關 : 細胞增生 (Proliferation) 、分化 (Differentiation )、存活(Survival)、自體凋零(Apoptosis)27,28。FGFRs 的活化主要是藉由跟細胞外的纖維母細胞生長因子(FGF, Fibroblast Growth Factor)結合而產生雙聚合 (Dimerization)並且在蛋白質上特定的酪胺酸產生自體磷酸化(Autophosphorylation)。藉由自體磷酸化刺激下游的蛋白質進而產生生物訊息的傳遞(Biological signaling pathway)。與 FGFR 相關的訊息傳導路徑中，主要可分為三個部分: Ras/MAPK、PI-3 Kinase、PLCg29,30。. 訊息傳導路徑的活化與細胞的種類及分化的階段有關，使的下游特定蛋白質標的活化 31。從 Figure16 中不難發現上述提到的三條訊息傳導路徑: 膜上的 FGFR1 區域先與 Frs-2 結合，結合後會使得 Grb2 活化並與 Sos 結合，之後身為 G-protein 的 Ras 會與 Sos 結合，最後活化 Raf/MAPK pathway。經研究指出，這條訊息路徑與 DNA 合成、增生以及分化有關 32；在 FGFR1 與 Frs-2 結合且活化 Grb2 後，Grb2 的 SH2 domain 會與 Gab1 結合，Gab1 使得 PI3-K 進一步的刺激 Akt 而導致細胞的存活. 50.

(53) 33. ；最後一部分是 FGFR1 的酪胺酸蛋白激酶區(TK2)上的 Tyr766 自體磷酸化並活化. PLC 蛋白質，PLC 的活化產生了 DAG 與 IP3 兩種化學物質，這兩種化學物質分別造成 PKC 的活化以節細胞內鈣離子的釋放 34。除了以上三種路徑之外，我們還可以發現 FGFR1 蛋白質與細胞增生的 Crk、與細胞骨架重新排列的 Src、肌動蛋白運動的 Nck/Pak/Rac 有關。. Figure 16: FGFR1 在生物細胞訊息傳導路徑所扮演的角色 35. 在最近的研究中，FGFRs 家族被發現與癌症有密不可分的關係。FGFRs 家族被發現在血管新生(Angiogenesis)、細胞型態改變(Changes in cell morphology)、細胞移動(Cell motility)及腫瘤增生 (Proliferation) 有關 36。特別是 FGFR1 的過度表現，會讓許多類型的癌症活躍。這些癌症包含腎、肝、肺、胰臟、大腸、皮膚、卵巢、前列腺以及血癌等等. 33,35,37. 。並且在最近的標靶藥物研發中，FGFR1 也被視為熱門的. 標的蛋白質。科學家們認為發展 FGFR1 的抑制劑為治療癌症上不可或缺的一環。. 51.

(54) 在結構上，纖維母細胞生長因子為一類型的膜蛋白，其結構 38 可分為 1. 接收細胞外訊息的三個免疫球蛋白區域(Ig-like domain)、2. 穿越細胞膜區(Transmembrane domain)、3.兩個免疫球蛋白間的區域(Acidic box between IgI and IgII)以及 4. 細胞內的與 ATP 結合的酪胺酸蛋白激酶區(Cellular tyrosine kinase domain)。由於酪胺酸蛋白質激酶區的活化情形與疾病最有關連，故在我們的研究中，僅討論此區中 ATP 結合的活化位置與化學小分子抑制劑的關聯。. Figure 17: FGFR1 酪胺酸蛋白質激酶區的 T 結構(pdb code:3JS239). Figure 18: 3JS2 蛋白質序列(綠: Beta Sheet；黃: Alpha Helix). 52.

(55) 1.3. 研究目標 (Research Objective) 在以下的研究中，我們使用分子嵌合模擬的方法針對 FGFR1 酪胺酸蛋白質激酶. 區，探討在不同模擬設定條件下搜尋出分子資料庫中具有活性分子的比率為何，並且最佳化設定條件，以用來提升虛擬篩選研究時找到具有生物活性化合物之機率。. 53.

(56) Chapter 2. 模擬方法 (Simulation Methods). 分子嵌合模擬 (Molecular Docking). 2.1. 在分子嵌合模擬中，我們使用較簡略的評分函數取代複雜度較高的力場模型來描述原子間作用力，因而可以快速的篩選大量的分子資料庫中具有潛力成為抑制劑的化學結構。關於我們所使用的軟體及評分函數，將在下面做更進一步的介紹。. 2.1.1. GOLD 與基因演算法用來做分子嵌合計算的軟體有很多，如Dock40、AutoDock41、Glide42、ICM43、. GOLD44等。在我們的研究中所使用的是GOLD軟體。GOLD（Genetic Optimisaton for Ligand Docking）是由英國Cambridge Crystallographic Data Centre所開發的一套專門用來做分子嵌合的軟體，此套軟體使用基因演算法計算蛋白質與藥物分子之間最佳化的結合構型，並在計算時會讓小分子上的化學鍵轉動並可在結合位置附近放入溶劑分子。這些動作是為了尋找小分子結合在蛋白質活化位置附近時所呈現的最低自由能，以便找到兩者之間結合時正確分子構型。. 在達爾文的演化論中提到『物競天擇，適者生存』，物種在不斷變遷或惡劣的環境中為了生存及適應環境，而不斷演化，產生適應力更好的下一代。生物在繁殖的過程中染色體會進行交配及突變來改變基因的組成，使得子代和親代及子代之間產生差異性。生物就是利用這種繁殖的機制，造成無數的變異發生，使得每一個個體有著不同的特性。在自然環境的考驗及彼此之間的競爭下，將不適生存者淘汰，而生存者繼續藉由繁殖。在這之中不僅將優良的基因延續，還有可能使下一代擁有更適應環境的基因。. 在基因演算法45中，我們稱每一個體稱為染色體。每一染色體的數值是由亂數產生，而每一世代的染色體所成的集合稱做群體（Population）。在每一世代中的每一個染色體互相競爭，較適合生存環境的則有較高的適性值（Fitness value），而有較高的適性值的染色體可以複製出較多的子代，然後從其中選擇配對來交配產生下一代，以產生適存度更高的下一代。再者，並加入突變的處理，以產生更多的變異。按：在基因演算法中評分函數（Scoring Function）稱為適性化函數（Fitness Function） 54.

(57) 如此一代一代的演化下去，將產生適性值較高的染色體，而該染色體即是我們找到的解。基因演算法的限制為在搜尋解答時不保證找到真正的最佳解(extreme minima)，而是找出局部最佳解(local minima)，但此局部最佳解是經過廣大解答空間演化到某種程度的可能最佳解。對於若干無法預知最佳解的狀況下，此演算法可以快速求得某種程度之滿意解。因此適合非線性等多變數複雜問題之求解。遺傳演算法的基本運算為選擇、交配和突變等機制，而這些機制之操作狀態通常是根據隨機值而改變的，因此即使在環境參數完全固定不變下，每一次運算求解所得之答案可能並不相同。基因演算法流程圖如 figure 19:. Figure 19: 基因演算法之流程圖 45. 在基因演算法中，Number of Operation 參數所指的是達到收斂前選擇、交配和突變三種操作的最大值，在我們的研究中也嘗試著更改不同的 Number of Operation 來驗證整體研究的收斂性。 55.

(58) 評分函數. 2.1.2. 在GOLD中的評分函數分為下列幾種：GoldScore46、ChemSore46、ASP46 (the Astex Statistical Potential)。在實驗室先前的討論中，發現ChemScore比較適合針對蛋白質激酶(Kinase)做研究，因此這邊我們所使用的評分函數皆為ChemScore47。評分函數大概可分為三個種類48: Force Field、Empirical、Knowledge-based。而 ChemScore屬於Emperical類型的評分函數。其函數形式如下:. ∆ +∆G. =∆. +∆. (∆r). (. )+∆. (. (∆α). )+∆. Eq. 9. 在函數中G的係數為未知且必須從線性回歸所取得。在這個總共可以分為四個部分: 1. Ghbond代表的是氫鍵的作用力在氫鍵作用力中分為acceptor、donor。所謂的acceptor指的是沒有與氫原子鍵結在一起的氮原子、氧原子末端的硫原子及離子態的鹵素原子；相對的donor為接有氫的氮、氧原子。g1、g2為係數，這兩個係數分別跟距離與角度有關，其仔細描述如下:. 1, if ∆r ≤ 0.25Å (∆r) = 1 − (∆r − 0.25)⁄0.4 , if 0.25Å < ∆r ≤ 0.65Å 0, if ∆r > 0.65Å 1, if ∆α ≤ 30° (∆α) = 1 − (∆α − 30)⁄50 , if 30° < ∆α ≤ 80° 0, if ∆α > 80°. Eq. 10. 2. Gmetal代表的是金屬與其他原子的作用在這裡我們希望的不只是能夠描述有機分子的評分函數，並且在許多蛋白質的活化位置也有著重要的金屬原子存在。因此在這裡，ChemScore 特別將金屬的作用力提出來分別描述其與非金屬原子的關係。之中 56.

(59) ∑. (. )所代表的是所有acceptor/donor與分子內金屬由距離所產生的係. 數大小。如Fig所示，R1為2.2Å、R2為2.6 Å。. 3. Glipo為親脂性官能基的作用親脂性官能基在ChemScore中的定義為:氯、溴、碘、硫、碳原子。之中∑. (. )所代表的是所有親脂性原子I、L間距離所產生的係數大小。如. Fig所示，R1為. +. +0.5、R2為(R1+3.0) Å。其中. 、. 為I、L兩原. 子的凡得瓦半徑。. Figure 20: Gmetal、Glipo 中 f(r)與 R1、R2 之關係 47. 4. Grot為可旋轉的鍵之間的作用最後一項代表的是辨別可旋轉的化學鍵。在這邊對於可旋轉鍵的定義為sp3-sp3或sp2-sp3種類的化學鍵。但這些化學鍵不包含末端的CH3、CF3、 NH2、NH3。只有在化學鍵兩邊都與receptor接觸時，才會被認為是不可旋轉的鍵；但若化學鍵有一邊可以旋轉而另一邊不能旋轉時，這個時候熵值的貢獻會顯得明顯。因此這邊我們特別把這一項列出來。下列的函數用來表示固定化學鍵時對擺動的貢獻:. = 1 + (1 − 1⁄. ). (. ( )+. ( ))/2. Eq. 11 在這邊 Ntot 指的是不可旋轉鍵的個數，Pnl(r)、P’nl(r)為非親脂性重原子在化學鍵兩端所佔的比例。 57.

(60) 2.1.3. 豐富指數與 DUD 資料庫. 豐富指數(Enrichment. Factor, EF)49 是用來評估虛擬篩選表現的一種指標。在豐. 富指數中，我們將分子分為兩大類:第一大類稱為 decoys。所謂的 decoy 指的是未知活性的化學分子；第二大類稱為 active ligands，這一類型的分子為已經驗證過對特定蛋白質有抑制活性的分子。在豐富指數的實驗中，我們會同時將 decoys 及 active ligands 放入篩選的名單中，然後設定不同的模擬條件(如胺基酸主鍊或側鍊的擺動、蛋白質氫鍵位置的限制)。若在篩選中，我們所設定的模擬條件可找到比例較高的 active ligands，則代表此種條件較適合用來作為進一步的虛擬篩選設定。豐富指數的定義如下:. × Eq. 12 在這裡 Hitssel 指的是利用電腦篩選出的 active ligands 數量；Hitstot 為我們放入的 active ligands 的數量；NCtot 為 active ligands 及 decoys 的總數量；NC 為某個特定百分比內所有分子的數量。藉由豐富指數，我們可以得知用 docking 方法所挑選出來的 active ligands 是否優於隨機所挑選的 active ligands 數量。若此豐富指數為 1，則代表使用 docking 方法挑選 active ligands 的效率與隨機挑選的效率相同。換句話說，若豐富指數為 10，則代表使用 docking 方法挑選 active ligands 的效率為隨機挑選的 10 倍。. 在做豐富指數的實驗時，我們所採用的 active ligands 與 deocys 為美國 UCSF 藥物化學系 shoichet 教授所製作的 DUD 資料庫(A Directory of Useful Decoys)50。在 DUD 資料庫中，總共提供了 95316 個 decoys 分子、2950 個 active ligands，分別針對了 40 個不同的蛋白質標靶。而這邊我們使用 FGFR1 蛋白質的 ligands / decoys 組合。在 DUD 中 FGFR1 Kinase 的資料庫中，active ligand 的個數為 120 個分子且 decoys 的數目為 4550 個分子。但是為了畫圖方便，在豐富指數部分的實驗我們所畫圖的方式為在特定比例下挑選出 active ligands 的數量。. 58.

(61) 2.2. 實驗設定 (Experimental Setting) 在我們的研究中，我們所使用的蛋白質結構為蛋白質資料庫中人類的纖維母細. 胞生長因子受體一。在我們的研究中，我們設定蛋白質條件如下:固定蛋白質、單獨擺動活化位置 10Å 內胺基酸側鍊(Side-chain)、限制氫鍵位置(在這邊我所限制的氫鍵位置為 Ala56439 主鍊上面醯氨基(-CONH-)的氧原子及上的氫原子)、同時擺動胺基酸側鍊及限制氫鍵位置。在擺動胺基酸測鍊的部分，我們並沒有讓選取的氨基酸完全自由擺動，我們所利用的為 S. C. Lovell 在 2000 年發表的側鍊擺動資料庫 51(Rotamer Library)。我們利用上述的條件做豐富指數實驗，分別評估在前 1%、5%時的豐富指數、相對於蛋白質固定不動(rigid)的豐富指數以及挑選出的數量。. 我們先利用 EasyModeller52,53 軟體，將原本 DUD 資料庫中 FGFR1 所採用的結晶結構 1AGW 做 homology modeling ，將完整補齊 1AGW 的序列另存為 1AGW-complete。並且使用 2010 年 Krishna P. Ravindranathan 等人發表在 J.Med. Chem. 上面的結晶結構 3JS2 中的 a、b chain 的兩個結晶結構。我們一開始利用這四個序列完全相同但三度空間位置不同的結晶結構搭配上一段所提到的設定條件以及 GOLD 程式中不同的 Number of Operation 做以下的研究。在研究的過程中，我們發現 3JS2-b chain 更適合用來做 docking 研究，因此研究的後半部是以 3JS2-b chain 為主體。在我們的研究中為了使整體問題簡化，所以在 docking 實驗中並沒有加入結晶結構中的水分子。. Figure 21: FGFR1 蛋白質序列比較:1AGW1 及 3JS2. 59.

(62) Chapter 3 3.1. 結果與討論 (Results and Discussion). 側鍊擺動在 Docking 的第一個部分，我們先將上述四個結晶結構針對 DUD 資料庫中的. active ligand/decoys 做測試，並且擺動結晶位置中心 10Å 內的殘基(Residues)上面的側鍊(Side-chain)，以觀察是否擺動側鍊時比較容易挑選出 DUD 資料庫中的 active ligand。Figure 22 是以 3JS2 b-chain 為例，在 Figure 22 中我們擺動了 3JS2 b-chain 結晶位置 10Å 內共 28 個側鍊，並且觀察在前 1%((120+1550)*0.01, 前 46 個分子中)、 5%((120+1550)*0.05, 前 233 個分子中)所挑選到 active ligand 的數量，在這邊我們所使用的 Number of Operation 為 GOLD 程式中的預設值 10% (3,000):. Figure 22: 擺動 3JS2 b-chain 中 10Å 內側鍊擺動與 active ligand 的挑選數量分佈圖 60.

(63) Figure 22 中可以發現，擺動大部分側鍊時，所挑選出來具有活性的分子數量並不會比固定蛋白質(Rigid)條件下增加太多，甚至有些特定的側鍊在擺動後所挑選的數量反而是降低的。但是可以注意的是在擺動 514 號胺基酸(Lysine)時，其增加的挑選效果明顯優於其他的側鍊，也就是說擺動這個側鍊時更能挑選出具有活性的分子。從 Table 17 中可以更明顯的看到這個效應。. Table 17: 擺動 3JS2 b-chain 中 10Å 內側鍊擺動與 active ligand 的挑選效率 Condition rigid 482 484 486 491 492 493 494 511 513 514 545 561 562 563 565 566 568 569 571 627 628 629 630 631 632 638 637 641. No. of Ligands 1% 5% 2 8 2 6 1 13 3 11 1 8 1 11 3 8 2 9 1 10 4 12 12 29 4 9 3 9 1 13 1 9 5 11 2 11 1 7 6 9 2 9 1 10 4 8 4 9 6 11 3 8 3 10 4 9 3 10 4 8. EF 1% 1.66 1.66 0.83 2.48 0.83 0.83 2.48 1.66 0.83 3.31 9.94 3.31 2.48 0.83 0.83 4.14 1.66 0.83 4.97 1.66 0.83 3.31 3.31 4.97 2.48 2.48 3.31 2.48 3.31. 61. 5% 1.33 1.00 2.16 1.83 1.33 1.83 1.33 1.50 1.66 2.00 4.82 1.50 1.50 2.16 1.50 1.83 1.83 1.16 1.50 1.50 1.66 1.33 1.50 1.83 1.33 1.66 1.50 1.66 1.33. EF Ratio to Rigid 1% 5% 1.00 1.00 1.00 0.75 0.50 1.63 1.50 1.38 0.50 1.00 0.50 1.38 1.50 1.00 1.00 1.13 0.50 1.25 2.00 1.50 6.00 3.63 2.00 1.13 1.50 1.13 0.50 1.63 0.50 1.13 2.50 1.38 1.00 1.38 0.50 0.88 3.00 1.13 1.00 1.13 0.50 1.25 2.00 1.00 2.00 1.13 3.00 1.38 1.50 1.00 1.50 1.25 2.00 1.13 1.50 1.25 2.00 1.00.

細胞外調節激酶、纖維母細胞生長受體激酶與其抑制劑結合自由能之計算： 熱力學積分分子動態及分子嵌合模擬研究

細胞外調節激酶、纖維母細胞生長受體激酶與其抑制劑結合自由能之計算：熱力學積分分子動態及分子嵌合模擬研究