穿膜蛋白質之結構預測

(1)

第一章緒論

1-1 研究動機與目的

穿膜蛋白 (transmembrane proteins，TM proteins) 在細胞與細胞間訊息與物質的傳遞過程裡扮演著相當重要的媒介，同時掌控著我們的感官系統。若其功

能失常將導致多種疾病，此家族的部份成員甚至是目前市面上所見，將近一半

藥物所鎖定以發揮作用的對象 (Fleming,2000)。雖然已經由實驗決定結構的蛋白質數量與日俱增，但其主要為水溶性蛋白。穿膜蛋白種類雖佔全部蛋白質的

20-30 % (Krogh,2001;Liu & Rost,2001)，但時至今日，根據 UC Irvine 的 Stephen White 實驗室的統計，具已知結構的穿膜蛋白尚不到 190 種 (http://blanco.biomol.uci.edu/Membrane_Proteins_xtal.html#Latest)。而其主要原因為：要同時使膜蛋白完好的呈現原始折疊結構，又要使其變成具高度純化的

結晶體，一直仍是X-ray 晶體繞射 (X-ray diffraction crystallography) 技術的一大門檻。

對於蛋白質摺疊模擬而言情況正好相反：水溶型蛋白因其具有多種自由

度，在電腦模擬中可能出現相當多可能的摺疊態；但就在真核生物的細胞膜內

發現的穿膜蛋白而言，當由胺基酸構成的多肽鏈 (polypeptide chain) 處於雙層磷脂 (phospholipid bilayer) 的非極性區域時，會將較具極性的骨幹 (backbone) 的部分藏起，而露出相對較多厭水性 (hydrophobic) 的旁鏈 (或稱 R 基， sidechain)，而呈現 α 螺旋 (α helix) 的二級結構。因為雙層膜性質上的限制而減少了自由度，並且此種蛋白質的三級結構便是由數根α 螺旋聚集而成 (另外尚有β sheet的二級結構組成的β barrel型態之穿膜蛋白，但只在原核生物的外膜中被發現)，因此相較於水溶性蛋白，穿膜蛋白降低了電腦模擬預測的複雜度。即使如此，穿膜蛋白的模擬至今仍然是一個很大的挑戰 (Bowie,2005)，已有很多的學者在穿膜蛋白的模擬上做過相當多的嘗試。Mitaku 和 Hirokawa 等人針對預測Bacteriorhodopsin 的原始結構發展了一個將穿膜螺旋視為連續柱體

(2)

的方法。雖然部份的α 螺旋位置被成功的預測出來，依然不是原子結構等級的預測結果 (Mitaku et al.,1995;Mitaku et al.,2000)；Goddard III 等人則發展了預測 G 蛋白連結型接收器 (G protein-coupled receptors) 結構的模擬方法。然而他們

的模型需建立在bovine rhodopsin 的電子密度圖上，且結構最佳化的方法也被

限制在以電子密度圖原始分布作方向角的旋轉上 (Goddard III et al.,2002)； Fleishman 和 Ben-Tal 等人特別針對殘基的親、厭水性，穿膜蛋白家族胺基酸序列裡最常出現的幾個殘基，特別容易出現在不同螺旋間接觸面上的特性等等，

發展出一套計分函數 (score function)，以此預測出原子等級的 rhodopsin 結構圖。然而其模型同樣需要建立於低溫電子顯微鏡 (cryo-EM) 的低解析度初始結構上，且模型自由度也被限制於螺旋的方向角轉動 (Fleishman et al.,2004)。在穿膜蛋白家族中的眾多成員之中，我們選擇了 retinal proteins 的 bacteriorhodopsin 和 halorhodopsin 來作為我們的研究目標。因其已有高解析度的X-ray 結構圖，可使我們方便檢驗模型的效度。並因其原始結構呈現多肽鏈依序環繞的構形，相較之下模擬複雜度較低，極適合做為穿膜蛋白模擬模型的研究起點。我們期望可以建立一個，只需要以胺基酸序列作為先決條件的蛋白質摺疊模型。以較為擬真的連續模型 (off-lattice model)，用 Cα原子為骨幹的真實螺旋結構為二級結構，計算以retinal proteins 為藍圖發展的系統自由能變化，來得到可與其原始結構作比擬的原子結構等級之預測結果。 1-2 研究架構我們的研究方法建立在三個理論之上：一為 Metropolis 等人於 1953 年發表的以蒙地卡羅法 (Monte Carlo method，或 MC method) 為基礎作增強的 Metropolis 演算法 (Metropolis et al.,1953)，並結合 Swendsen 與 Wang 在 1986 年提出的replica exchange (以下簡稱 RE) 演算法 (Swendsen & Wang,1986)。二為Popot 和 Engelman 所提出的 two-stage 模型 (Popot & Engelman,1990)。最後則為Anfinsen’s dogma (Anfinsen,1973)。我們根據以上理論將研究分為以下幾

(3)

項架構流程：

(1) 將目標蛋白質上各穿膜 α 螺旋的胺基酸序列放入 AMBER9 force field 2003 執行能量最小化 (minimization of energy，以下簡稱 ME)，來得到所需的穿膜二級結構。

(2) 以(1)所得到的結構建立以 Cα為骨幹之連續模型，並將一系列的非共價

性作用反應置入模型中，以計算系統自由能。

(3) 以 Metropolis 演算法內的採樣法則來決定系統的演化方向。 (4) 若干 MC steps 之後執行 RE 來加速尋找具較低能量構形的過程。 (5) 重複(3)與(4)直到找到具整體最小能量的構形 (global energy minimum

structure，以下簡稱 Emin 結構)，此為模型之下最為穩定的結構。

(6) 將所得之 Emin 結構與已知的 X-ray 晶體繞射結構做分析比較。

1-3 目標蛋白質

Retinal proteins (以下簡稱 RPs) 為一與 retinal 分子互相鏈結，且具光子接受器功能的蛋白質，可在動物的視網膜上或是古細菌的細胞膜上被發現。 Retinal 是維他命 A 三種形式中的一種，為一多烯基發色團 (polyene chromophore)。其尚未與蛋白質產生鍵結前的分子式為C H O 。當它欲與20 28 RPs 的第七根 α 螺旋上的殘基 lysine 共價鍵結時，其上的碳原子便會損失一個氧原子， lysine R 基上的氮原子則損失兩個氫原子，而形成C N= 雙鍵鍵結 (此即所謂Schiff base)。一般而言典型的 RPs 存在於動物的視網膜上，其中一個有名的例子像是 rhodopsin，它同時也是 G 蛋白連結型接收器的成員之一。而我們研究的範圍則是存在於原核生物內的RPs。兩種類型的 RPs 都是七根以六段鏈圈 (loops，通常由數個親水性殘基串連而成) 連接而起的穿膜 α 螺旋鏈，相聚而成如桶狀般的結構。立體構造端見示意圖 Fig.1-1。多肽鏈的 N 端伸出膜外暴露在細胞外的環境，C 端則是在細胞質內，除此之外，兩類型的 RPs 在胺基酸序列相似度極低的證據下，並沒有任何其他可說明他們在演化上的相關性

(4)

之佐證。存在於原核生物中的 RPs 包括了 bacteriorhodopsin、halorhodopsin、 proteorhodopsin 和 sensory rhodopsin，而前兩項便是在本次我們 RPs 的研究目錄內。

Fig. 1-1 RPs 在雙層脂質膜中的結構示意圖。圖中為嗜鹽性古細菌 Halobacterium salinarum 內的四種 RPs。由左至右分別是 bacteriorhodopsin、halorhodopsin、sensory rhodopsin I 和 sensory rhodopsin II。

本圖摘錄自 (http://www.biochem.mpg.de/en/rd/oesterhelt/web_page_list/Topic_retinal_proteins_ Hasal/index.html) 1-3-1 Bacteriorhodopsin (菌視紫紅質，以下簡稱 BRD) 當嗜鹽性古細菌Halobacterium salinarum 處於低氧的環境時，它的細胞膜上會長出一塊塊的紫斑區域，即所謂purple membrane。其幾乎可佔古細菌表面的50%，而 BRD 便是在此區域被發現。它也是 purple membrane 上唯一的蛋白質，功能是作為一個光驅動的質子幫浦 (light-driven proton pump)。在正常的情形下吸收光子後，其內的 retinal 分子在產生異構化的過程，會使得質子從細胞膜在細胞質內的一端被傳遞到細胞膜外的另外一端，來對抗原本的濃度梯度所造成的質子傳輸方向性。而在過程裡BRD 內

部截面的residues 種類決定了質子的傳送路徑 (Kimura et al.,1997)。簡而言之便是它可將綠光 (波長介於 500-650 nm，但 568nm 可達最大吸收功率) 的能量轉換為電化學性的質子梯度，把質子傳送到細胞膜外。這些質子會參

(5)

時氧氣所扮演的角色。另外BRD 同時也利用產生膜電位改變的方式，涉及趨光性的運作 (Bibikov et al.,1993;Armitage et al.,1996)。

BRD 位於 purple membrane 中，天生適宜處在低氧濃度與高光密度的環境之中。這樣天然的二維晶體結構提供給了世人一個得以窺視穿膜蛋白的機會，也使它成為了研究目標的焦點之一。從1990 年 Henderson 等人製造出了第一個原子等級的BRD 結構圖 (Henderson et al.,1990) 後開始，解析度越來越高的結構圖紛紛問世。在Luecke 於 2000 年所發表的回顧論文裡， PDB 內便已刊登了 20 個 BRD 的蛋白質原子級結構座標 (Luecke,2000)，其實驗方法同時包含了電子繞射與 X-ray 晶體學。近年關於 BRD 的研究重心，一直在於上述的光驅動性質子幫浦運作背後的真正力學機制為何。其研究幕後功臣便是 PDB 內可詳細分辨 α 螺旋鏈上，各個 R 基、水分子與 retinal 如何相互作用的高解析度結構圖。1997 至 1999 的兩年時間， BRD 結構圖解析度便從2.9Å 進步到 1.55Å。其中包含了 1.55 Å 的 1C3W (Luecke et al.,1999)、1.9 Å 的 1QHJ (Pebay-Peyroula et al.,1999)、2.3Å 的 1BRX (Luecke et al.,1998)、2.5Å 的 1AP9 (Pebay-Peyroula et al.,1997)和 2.9Å 的 1BRR (Essen et al.,1998)，以上的英文數字串皆為 PDB code。在我們的模型裡被供作模擬結果比較之用的，為較近年的實驗結果1PY6 (Bowie et al.,2004)，以下為

其PDB 上之基本資料與 3D 結構圖。全名：以bicell 法結晶而成的 Bacteriorhodopsin。發表日期：2003/12/16 (最後修正 2009/02/24)。實驗方法：X-ray 繞射。解析度：1.8 Å。蛋白質分類：膜蛋白。分子來源：Halobacterium Salinarum。分子類型：多肽鏈 (共 249 個殘基)。

(6)

配體種類：retinal (C20H28O)。

Fig. 1-2 BRD 之原始 3D 結構圖。採 PDB code = 1PY6 之結構數據，配合軟體 MOLMOL 繪出。在本論文第四章的 BRD 結果討論中，皆以此結構作為比較基

礎。以細胞質內往細胞外之方向視角。左端最為突出的α 螺旋為 H1，其餘六根

以順時針的方向排列環繞retinal (圖中沒有繪出)。

1-3-2 Halorhodopsin (嗜鹽視紫紅質，以下簡稱 HRD)

同樣來自於古細菌Halobacterium salinarum 的 HRD 是另一種光驅動性離子幫浦 (light-driven ion pump)。它一樣可將綠光轉換成能量對抗膜的電

位能以及濃度梯度來傳輸氯離子進入細胞內，但和BRD 不一樣的是具有最佳吸收功率的光波長為 578nm。除了作為氯離子的幫浦，它同樣也可傳輸其他鹵化物或硝酸鹽等等。此種陰離子幫浦存在的目的，主要是因為此種嗜鹽性的古細菌，當身處於與本身細胞質的鹽分濃度等滲透壓的濃鹽水環境中時，便會開始急速生長。生長過程中如何維持自己本身的滲透壓，對細胞生存而言為一必要條件。因此他們必須不斷的自環境補充氯化鉀分

(7)

子，而這種光驅動的陰離子幫浦便可以替它們省下不少新陳代謝所需的能量。就目前的研究而言，HRD 是唯一個被發現的靠光驅動之陰離子幫浦 (Essen et al.,2000)。HRD 的性質與功能也被與 BRD 作大量的比較以釐清他們之間的相異之處，他們之間最大的差異，在於同樣身為光驅動的離子幫浦，他們所傳輸的離子具相反電性並且傳往相反的方向。就如同其他來自於Halobacterium salinarum 的 RPs 一般，HRD 以七根被六段鏈圈串接而起的穿膜 α 螺旋鏈，折疊而成一桶狀結構。H1-H7 緊密的環繞著retinal 分子，其以 Schiff base 與 H7 上唯一的 lysine (Lys242) 作共價鍵結。且retinal 分子把由 H2、H3、H6 和 H7 所圍成的通道分隔成面細胞內與面細胞外兩個區域。和BRD 之間雖有上述的功能性差異，在結構上則具有高度的相似性，在胺基酸序列的相似度上也達 31%。特別是 H3-H7 的這五段 α 螺旋共 122 個 Cα原子的結構座標，在均方根偏移 (root mean square deviation，以下簡稱 RMSD) 的表現上，兩者只有 0.74Å 的差異。若將兩個蛋白質結構以H3-H7 的區塊為準軸疊合起來，則會發現在 H1 和 H2 兩者的RMSD 分別是 2.1 和 1.3Å。將以上數據匯整而出的，便是成功地在

立方晶型的油脂面 (cubic lipidic phase) 中，將 HRD 長成晶體結構，並用 X-ray 繞射得到解析度高達 1.8Å 結構圖 (PDB code 1E12) 的 Oesterhelt 等人 (Essen et al.,2000)。2006 年 Oesterhelt 的研究團隊，又以 HRD 的突變種 T203V 作觀測得到其穩定態結構圖 (PDB code 2JAF)，解析度達 1.7Å (Tittor et al.,2007)。在我們的研究中選擇了 1E12 作為結構比較之依據，以下為其資料與結構圖。全名：Halorhodopsin，光驅動性氯離子幫浦。發表日期：2000/06/02 (最後修正 2009/02/24)。實驗方法：X-ray 繞射。解析度：1.8 Å。

(8)

蛋白質分類：離子幫浦。

分子來源：Halobacterium Salinarum。分子類型：多肽鏈 (共 253 個殘基)。配體種類：retinal (C20H28O)。

Fig. 1-3 HRD 之原始 3D 結構圖。採 PDB code = 1E12 之結構數據，配合軟體 MOLMOL 繪出。在本論文第四章的 HRD 結果討論中，皆以此結構作為比較基

礎。以細胞質內往細胞外之方向視角。左端最為突出的α 螺旋為 H1，其餘六根

(9)

第二章演算方法

2-1 蒙地卡羅法 (MC method) 2-1-1 歷史簡介蒙地卡羅的名字，主要是借自於摩納哥 (Monaco) 境內一個以賭場而聞名世界的城市。賭場內常出現的一種輪盤遊戲 (Roulette)，其以產生亂數且不停重複的遊戲特性，與蒙地卡羅法的核心概念正好不謀而合。以亂數為基礎的計算與實驗方法的起頭，至少可追溯至19 世紀末期。初期運用此方法最著名的例子或許可說是 1930 年代的 Fermi。他用了亂數方法來計算那時剛被發現的中子特性。但他並沒有發表任何與他所使用的統計方法相關的論文，因此也沒有給予命名。蒙地卡羅法命運轉折的里程碑則是 1945 年。在那個極具科學爆炸性的年度裡，發生了兩件科學史上的大事：一是在新墨西哥州的Alamogordo 舉行的第一次原子彈試爆成功；二則是第一台

電子計算機 (Electronic Numerical Integrator & Computer，以下簡稱 ENIAC) 的誕生 (區別於 Zuse 在 1941 年架設的 Z3，其是利用電磁而非電子)。兩件乍看之下並無相關的事件，卻在兩派人馬的晤面下碰出了火花。出席的科

學家中包括了Fermi、von Neumann、Ulam 和 Metropolis。雖然以亂數採樣技術為主幹的原子彈研發計畫獲得了空前的成功，但此項統計技術隨著計

算式的日益冗長而似乎走到了人力無法負荷的絕境。當上述的物理學家們

看到了那條被ENIAC 誕生所帶來的無限可能性照耀出的路時，那條路無疑

地便是通往蒙地卡羅法。二年之後，Ulam 和 von Neumann 提出了以 ENIAC 結合統計方法，解決中子在可裂變物質中擴散問題的方案。並在隔年得以

實際測試。Ulam 和 Metropolis 在 von Neumann 和 Fermi 的協助下，解出了設計問題下的薛丁格方程本徵值，也證明了此統計方法的有效度。因此在

1949 年他們正式以蒙地卡羅法為演算法的名字來發表論文 (Metropolis & Ulam,1949;Metropolis,1987)。在那之後，隨著電腦計算速度的飛奔馳騁，蒙

(10)

地卡羅法被廣泛的運用在物理科學的範疇內，其中又包含計算物理、物理化學與統計物理等領域。另外又在通訊、視覺設計以及經濟金融等的領域中扮演重要角色。 2-1-2 原理內容要了解蒙地卡羅法的基本精神我們可由一“計算 π＂的實驗來著手。假想牆上有一正方形區塊，其內有一半徑＝正方形邊長的1 4圓，如Fig. 2-1。有一玩家連續地將飛鏢射往正方形區域內。如果玩家可以限制射點完全在方形內，但不特別瞄準某處。那麼我們可以藉由統計所有射點在灰色區域內的飛鏢數與總射擊飛鏢數，來完成 Eq.2-1 的計算即可得到 π 值。 1 4 = = (2-1) 4 π 圓的面積灰色區域內的飛鏢數　　。黑色方框內的飛鏢數正方形的面積 Fig. 2-1 實驗標靶圖。但要得到準確度在可接受範圍內的π 值，總飛鏢數可能最起碼要超過 1000 以上。也因此電腦的發展才會和蒙地卡羅法的命運如此息息相關。因為電腦可扮演的便是飛鏢手的角色：一個可以不會疲勞，且永遠將飛鏢丟擊在設定範圍內的玩家。任意丟擲的飛鏢即所謂亂數，所做的計算便是機率統計。凡是以此概念所做的實驗都可稱作蒙地卡羅實驗。在1953 年 Metropolis 等人又發表了一演算法，其特別針對以電腦來運算涉及分子間作用能量的物質特性研究。此修正後的蒙地卡羅法，將依照波茲曼分布 (Boltzmann distribution) 下的機率因子，來選擇下一步的新構形以取代原本完全隨機的選法 (Metropolis et al.,1953)。後人將此演算法稱為 Metropolis 演算法，或直接仍以蒙地卡羅法稱之。更具體地來說，我們

(11)

所研究的系統將不會是孤立於環境中的，相反地它將與環境交換能量，而其交換指標將會是環境溫度 T。如果溫度越高，環境將會更願意將能量分給系統；但溫度越低時，環境便會讓系統處於一個低能量的狀態。我們可用波茲曼因子 (_e-E T_{) 來描述當溫度為 T 時 (此處使用自然單位制，即波} 茲曼常數 = 1k )，我們可發現系統能量為 E 之機率關係將是： -E T _(2-2) P∝e 　　　　　　　　　　　　　。接著它會產生一馬可夫鏈 (Markov chain)，即新的可能構形將只和目前的初始構形相關，和更早之前的構形則沒有相關性，且新的構形也將滿足於 Eq.2-2。如果在這個過渡期裡，未來的可能構形其擁有的系統能量為E ，_j 目前構形為E 。在溫度為 T 的環境下，兩者發生機率比： _i ( - ) -( ) (2-3) ( ) j j i i E T E E T E T P j e e P i = e = 　　　　　　　　。這個演算法特殊簡便之處便在於藉由 Eq.2-3，我們不需要真的知道隱藏在 Eq.2-2 中的比例常數究竟為何。Metropolis 等人 (1953) 提供了以下的步驟，來重複進行同時滿足精細平衡條件與Eq.2-3 的構形篩選動作。 (1) 自初始狀態 i 開始，其系統能量為E 。 _i (2) 從系統中隨機地選擇一個或一群粒子，並隨機地移動被選定的粒子。 (3) 計算改變後的狀態 j 所具有的能量E 。 _j (4) 計算這兩個處於過渡期前後的狀態能量差 ΔE = E_j

−

E_i。 (5) 若 ΔE < 0，便將狀態 j 作為新的初始構形，因為其具有更低的系 統能量 (根據 Eq. 2-2，此狀態應為相對之下系統較想成為的狀態之一)。 (6) 若 ΔE > 0，則以機率為 _e-ΔE T _{的情況下來判斷是否將狀態} _{j 設}

(12)

為新構形。也就是我們將取一隨機數γ ，且1> > 。如果γ 0 _e E T γ _≤ Δ _，則狀態 _{j 將取代狀態}_{i 成為下一任的初始狀態。反} 之，系統將回到狀態 i 來進行下一次的馬可夫鏈。從機率函數的形 式我們可以了解到，如果環境溫度越高，系統將越不在乎是否走往 “錯誤＂的方向。但如果溫度極低時，系統將被強制一貫地演化至局部能量最低的構形。我們可用一簡單的程式，來概括由狀態 i 演化至狀態 j 的接受機率為：

(

-

)

( ) min 1 , E T (2-4) P i→ j = e Δ 　，　Eq.2-4 又被稱為 Metropolis 準則 (Metropolis criterion)。另外亦可用 Fig. 2-2 來表示此演算流程。

(13)

Fig. 2-2 Metropolis 演算法之流程示意圖。初始態 i 初始能量E _i 新狀態 j 新系統能量Ej 計算新舊狀態間能量差 j i E E E Δ = − 0

E

Δ

>

Δ

E

<

0 取一亂數γ 1> > γ 0 _e E T γ _≤ Δ 狀態 j 取代狀態 i 成為新的初始態 _e E T γ _> Δ 狀態 i 不被取代 計算所需的物理量

(14)

2-2 Replica Exchange Method (或稱 Parallel Tempering) 2-2-1 歷史簡介自蒙地卡羅法 (以下此名皆指已包含 Metropolis 演算法部分) 問世以來，已成功地在各個領域裡解決了許多的問題，但他的成功卻不是全面性的。在諸如蛋白質摺疊或是玻璃質材料等的複雜性系統中，都遇到了無法有效率的逃脫出局部位能井的障礙。根據Eq. 2-4 可知，在越低溫的環境越容易發生此現象。這將使模擬的重點取樣只侷限於一小部分的結構空間，

而非全面。一個解決方法是結合replica exchange (以下簡稱 RE) 的模擬。 RE 最早是出現在 Swendsen 和 Wang 所發表的自旋玻璃模擬一文中 (Swendsen & Wang,1986)。在此文獻裡，其自旋玻璃的亂序系統，在一系列

不同溫度下的 replica 內進行模擬。然而在進行鄰近的 RE 時，只有部分的

構形資訊被交換。爾後較為我們所熟悉的，建立於全面構形資訊交換下的

RE 法是由 Geyer 於 1991 年提出 (Geyer,1992)。一開始這個新方法的應用

被侷限於統計物理的範疇內。後來在Hannsmann 運用此 RE 結合蒙地卡羅

(以下簡稱 REMC) 法的演算法於生物分子系統的模擬 (Hannsmann,1997)； Falcioni 與 Deem 也利用此法來幫助決定 X-ray 粉末繞射所得之圖譜結構 (Falcioni & Deem,1998)； Sugita 和 Okamoto 則是將 RE 與分子動力學結合，運用在蛋白質摺疊的模擬上 (Sugita & Okamoto,1999)。在此之後，RE method 的運用多元性便快速成長，其觸角也同時跨越了各領域，如物理、化學、生物、電機以及材料科學等等。然而，要進行RE 的模擬演算便意味著需要用上比單獨模擬時更多倍的電腦設備。這也是為何最初此方法的接受度不高的原因。但是近十幾年來，電腦硬體的淘汰更新率日新月異。且RE 法也被證實，其運算設備的需求度較高，反映到只需要花更少的運算時間，也促使其接受度大增。在可讓各

replica 獨立運算的 CPU 叢集 (cluster) 技術日漸提升以來，也讓 RE 法的效率表現較往日更為卓越。

(15)

2-2-2 原理內容 RE 的概括想法可由 Fig. 2-3 來表示。大圓表示一個二維的相空間。當模擬系統處於低溫時，則模擬產生的馬可夫鏈在此相空間裡所走過的“路徑＂，便容易像被侷限在圖中的灰色區塊一般，如同走到了相空間內的一小山谷內而無法輕易脫離該處。但當溫度高時，模擬系統便容易像充斥於圖中的箭頭一般，遊走於相空間內的每個角落。同時也有可能“路過＂所謂的灰色區域，這時也就有了交換的可能。而當交換發生時，低溫的系統便可趁高溫之勢，跳出當時的灰色區塊而到達相空間的另一個角落。根據 Metropolis 準則，低溫的模擬系統對於尋找最低能量，相對於高溫的系統將較具有效率，因為前者比較不會走往錯誤的方向。但也因為如此，找到的往往只是灰色區塊內的局部最低能量 (小山谷)，而被卡在其中花掉許多無謂的演化時間。因此在高溫系統的協助下，可幫助低溫系統確定究竟哪一個局部最低能量，才是真正的整體最低能量。一般而言每個replica 具有不同的溫度。誠如上面所述，我們必須讓不同的replica 之間構成交換的條件，才能真正的完成任務。此條件便是建立於不同溫度間其構形能量分佈曲線的重疊部分，如Fig. 2-4 所示。以一般

Fig. 2-3 二維的相空間示意圖。本圖節錄自 (Earl & Deem,2005)，但需注意在本章節我們對此圖的譬喻角度與原文有觀點上的不同。

(16)

E P(E) T4 T3 T2 T1 Fig. 2-4 四個不同溫度下 replicas 的高斯能量對機率分布圖。(T1 <T2 <T3 <T4) 電腦模擬常接觸的典型高斯能量分布模型來做為範例，當 1 2 3 4 T <T <T < ，如兩溫度間的 TT Δ 不要過大，其能量分布曲線便會產生重疊。也就是上頁所比喻的高溫可能“路過＂低溫的機會。重疊的部份越大，則代表其可能接受交換的機率越高 (Hansmann,1997)。假設我們在此 REMC 法下的模擬中，共有 M 個具有不同溫度的獨立 replica ，其中 1 2 M 1 M T <T <... T< − <T 。Ti 為第 i 個 replica 的溫度，E 為其系統能量，i

{

C C1, 2,....,CM

}

則為所有 replica 在相同時間點下各個的初始態構形集合。依循 Metropolis 演算法的基本核心，模擬流程將需要建立一馬可夫鏈。但在此改良過的模型裡，馬可夫鏈將需要從兩個步驟來看待： (1) replica 內部的構形演化。replica 間在做交換之前各是獨立運算的個體，因此內部所進行的新舊構形之間的交換率只跟自己的系統能量變化有關，且接受與否也依然遵守Metropolis 準則 (Eq. 2-4)。 (2) 第 i 個 replica 與第 j 個 replica 之間的構形交換。因為一般而言可 能進行的交換都發生於鄰近的溫度之間，故 j= + 。而這裡的交換i 1

(17)

意味著兩個 replica 構形資訊將全部被改變，即 old j new i C =C ； old i new j C =C 。並因為之間同樣涉及粒子間的構形變化，其間的交換機率同樣依循 Metropolis 準則。相異處則為須同時顧慮到兩群的粒子系統，故此其交換率為：

(

( ) ( )

)

(

( ) ( )

)

( ) ( ) ( ) ( ) (1 1 )( ( ) ( )) ( ) min 1, min 1, min[1, old new j i i i i j j j j i i j i i j j j i j i E C T E C T E C T E C T E C T E C T E C T E C T T T E C E C P C C e e e e − − − − − − + + − − ⎡ ⎤ → = ⎣ ⎦ = ⎡_⎣ ⎤_⎦ =

⋅

(

)

] = min 1,_eΔ Δβ E (2-5)_，其中 1 1 j i T T β Δ = − ，且 Δ =E Ej −Ei。從Eq. 2-5 可知道，其實並不需要特定限制 RE 只會發生在鄰近的溫度之間。因為這樣的選擇 1 1 1 i i T T β + Δ = − ，自然會是最佳的，其交換接受率將隨著Δ 的升β 高以指數關係下降。由以上兩部分我們可以勾勒出REMC 法的演算輪廓，Fig. 2-5 為其示意圖。首先讓每個replica 進行內部的蒙地卡羅演算。經過某特定量的 MC step 後，再讓 replica 間透過 Eq. 2-5 來決定彼此是否交換構形。而在這之前的 MC step 數對於此新型演算法的效率有一定的影響。因為其所代表的演算時間，需足夠讓低溫的 replica 找到當時構形的局部最低能量狀態，再與高溫的replica 互換構形，才有幫助演化的實質意義。但其演算時間又不需過長，以免在找到局部最低之後又浪費了過多無謂的等待時間，而拖慢整體演算效率。除此之外尚有決定鄰近 replica 的溫度差等的變因，都在決定整體演算效率上有重大的影響性，我們將在下一小節中另行討論。為了增進演算速度，我們使用了叢集的主從架構概念來進行平行運算，使用的工作站是由1.8GHz AMD Opteron™ 的四核心處理器串連而起。

(18)

運算過程將分為一個主 (Master) 與數個從 (Slaves) 兩個部份，且不論主或

從皆各佔一 CPU 硬體資源。示意圖端見 Fig.2-6。每個 Slave 即代表一個

replica，他們將在各自的溫度下進行蒙地卡羅運算。Master 則負責彙整各 replica 的工作，並計算每個 replica 是否與前後鄰居交換構形，但前後只能

隨機擇一。再將新的構形與能量分配給對應的 replica，以待下次的交換工

作。

Fig. 2-5 REMC 模擬法示意圖。每條直線皆代表一 replica 內構形變化，以 • 串起的實線是以 MC step (MCS) 為單位的時間軸，◎則為規律地穿插其中的RE step。 Fig.2-6 平行運算之主從架構示意圖。 Master 網路 Slave (Replica 1) Slave (Replica 2) ‧‧‧‧‧ Slave (Replica M)

(19)

2-2-3 重要變因

可使 REMC 法效率最佳化的變因大致包括以下三項：(1) 兩次 RE 之

間的 MC step 數目；(2) 模型中 replica 的總數；(3) 鄰近 replica 的溫度差 (Earl & Deem,2005)。

第一項變因甚少在文獻中被討論或通式化，原因可能為其主要和本身模擬所使用的模型粒子數多寡、能量分布曲線的範圍和低溫replica 內尋找到局部最低能量所需的運算時間有關，甚難將其全部納為物理量來討論。因此在我們的模型中，我們藉由模擬過程中 REMC 的效能表現來決定大小。例如：是否常發生高低溫replica 交換構形後，低溫 replica 尚未找到具較低能量的結構之前，構形又再度交換回到高溫處。或是運算時間較平均而言被明顯拉長 (高溫與低溫 replica 間所需運算時間相差過大，使低溫 replica 需在等待交換前閒置過久)。前者需拉長 MC step 數，後者則是縮減。第二項和第三項變因則可綜合討論之。以上節出現的 1 2 M 1 M T <T <... T< − <T 來說，一般而言，T 將會是系統模擬所感興趣的溫1 度，在T 出現最低能量構形的機率應最高。1 T 則並沒有所謂上限，但其至M

少必須讓其他所有replica 都不致困在局部最低能量內。Kofke 與 Kone (Kone & Kofke,2005)建議避免以固定的T 與1 T (最低與最高溫)，來找尋能使模擬M 效能最佳化的replicas 數的觀點。選擇在T 固定的情形下，考量預算中的硬1 體資源與系統粒子數，先來決定所需replica 的數量。其中若模擬的粒子數 N 增加，因平均的總能量大約會和 N 成正比，但能量分布的寬度大約只以 N 比例增加，使得能量分布曲線更加窄化。這將讓不同溫度間的能量分布重疊變少，使所需 replica 數也將大致隨 N 的比例增加 (Kofke,2002;Hansmann,1997)。再選擇能使每個相鄰的 replica 都可具有相等交換機率_{P C}( old →_Cnew)的溫度分布 (讓每個 replica 在每個溫度下都可幾

(20)

乎花上相同的模擬時間，也就是在溫度空間內構成自由的 “random

walk＂)。他們發現當在定容比熱(C )在所有模擬的溫度範圍內都設為定值v

的條件下，此時REMC 的交換機率_{P C}( old →_Cnew)將只與其溫度比例(

i j

T T )

有關。在記錄下每個replica 實際在所有溫度中“random walk＂所做的平均平方位移(MSD)對平均的交換機率P(_Cold →_Cnew)_{做圖後，發現無論模擬所}

設定的C 值為何，其 MSD 最大值皆發生在v P(Cold →Cnew)= 23%處。而因

曲線斜率與C 無關，故判定只要鄰近 replica 間不要有劇烈的v C 變化，即使v

整個溫度範圍的C 變化極大，上述結果依然可以通用 (Kofke,2002)。 v

然而，於實際的模擬中，可使P C( old →Cnew)固定的replica 溫度分佈

難以估算，因此多處於一機率範圍內的形式。文獻上交換機率的採用範圍

也自 20% − %30 (de Pablo et al.,2005; Mitsutake et al.,2004) ，到

40% − %70 (Kokubo & Okamoto,2009;Mitsutake & Okamoto,2004) 不等。溫度也自呈指數分布到近等比數列分布，即T Ti j ≅ 常數。

(21)

第三章模型

3-1 蛋白質二級結構的模擬

1990 年 Popot 和 Engelman 提出了膜蛋白摺疊的兩階段模型 (two-stage model) (Popot & Engelman,1990)。在這個模型裡，形成膜蛋白穩定結構的過程可分為兩個互相獨立的熱力學階段：(1) 各個穩定且不互相影響的二級結構將在雙

層膜中形成。這裡所指的二級結構絕大多數為α 螺旋鏈。(2) 螺旋之間開始互相

作用以形成更緊密及穩定的三級結構。此理論提出後被廣泛的接受並被大量的運

用在膜蛋白的模擬上 (Bowie,1997 ; MacKenzie & Fleming,2008)。在此一較為簡

化的概念下，我們在進行蛋白質摺疊模擬之前，先行建構模型中所需的初始α 螺

旋結構。有別於我們先前所使用的coarse-grained model (Chen et al.,2008;Chen & Chen,2009)，我們嘗試將序列相依 (sequence-dependence) 與真實螺旋導入我們的模型中 (而非先前的序列分組及棍狀結構)。因此我們擷取 PDB 中目標蛋白質上各穿膜α 螺旋的胺基酸序列，並使用 poly-ALA 螺旋 (φ= - 60°，ψ= - 40°) 作為初始結構，放入AMBER9 以 force field 2003 執行能量最小化 (ME) 以得到在所需胺基酸序列下擁有最低能量的α 螺旋結構 (Chen & Chen,2009)。因膜蛋白的

二級結構形成過程應在脂質膜中完成，我們將介電常數 (κ) 設為 2.0-3.0 之間

(Tsong,1990)。視殘基數目多寡調整最速下降法 (steepest descent method) 與共軛梯度法 (conjugate gradient method) 的步數，但兩者比例維持於 1：2 (Baginski et al.,2002) ，以得到最佳的二級結構。 Fig. 3-1 為其中一個範例結果 (Bacteriorhodopsin 的第一根穿膜 α 螺旋或 H1 ，其胺基酸序列為： EWIWLALGTALMGLGTLYFLVKGM)。隨後我們將所得結構中的所有 Cα座標取出，亦即圖中實心球體的球心部份，以此架構出骨幹的螺旋結構。並在第二階段的模擬中保持同樣的相對位置，亦即僅以α 螺旋為單位做移轉，螺旋內部構造維持不變。

(22)

Fig. 3-1 以 BRD 的 H1 序列模擬所得之 α 螺旋結構圖。實心球體代表 Cα原子，

其中深灰色部分為相對於淺灰色部分，具有較大或較鬆散的 R 基結構之殘基種類

(包含 Lys、Tyr、Trp、Phe、Ile、Leu 和 Met)。多肽鏈的 N 基位於圖右，紅色直線

依序串起相鄰的Cα，灰線部分則為R 基鍵結處。

3-2 以 Cα 原子為骨幹之連續模型 ( Cα - backbone off-lattice model )

依照 two-stage 模型，我們將第二階段的初始狀態設為，所有在第一階段中得到的穿膜 α 螺旋隨機的散布在雙層膜中，且初始的螺旋軸 (helical axis) 皆垂直的豎立在膜中，即平行於Ｚ軸的方向。順序相連的前後螺旋間頭尾的 Cα距離 (r ) 設有鏈圈長度限制，使其無法超過實際最大鏈長 (正比於螺旋間鏈圈所擁有_ht 的殘基數目，而兩相連的殘基距離約3.8 Å)。只要r 尚在鏈圈長度限制之內，αht 螺旋皆可自由活動。換言之，即此位能函數為一階梯函數。如r 小於相對應鏈長，_ht ( )=0rht Ρ ；反之則Ρ( )=rht ∞ ，因此也不可能通過蒙地卡羅法的機率函數考驗而成為新的構形。而如Fig. 3-2 所示，穿膜 α 螺旋可活動的空間分為三個區域：水－

雙層膜－水。一般而言膜厚約 30 Å (White & Wimley,1999;Popot & Engelman,2000)，但膜是具有流動性且容易隨著鑲嵌其中而互相作用的膜蛋白來改變其厚度的 (Konings et al.,1991;Popot & Engelman,2000)，故在我們的模型中

(23)

Fig. 3-2 穿膜螺旋於模型中的三個自由度示意圖。傾斜角(θ )、方位角(φ)及自轉角 (Ω )。Ｚ軸垂直於水平分布的雙層膜。圖中實心球體代表 Cα原子的位置，藍色球位於想像柱體的前方，灰色部份則位於背面。(螺旋與雙層膜未按比例繪出) 每一次的MC Step 裡將有一根 α 螺旋被隨機選取，並在滿足 ( )=0Ρ rht 的情況下作自由的移轉，其中包括：質量中心在ＸＹＺ三維空間的移動，並可改變螺旋本身的φ角與θ 角，另外因為螺旋亦可視為一表面具有可分辨性的圓柱體，所以亦具繞著自己的轉軸做自旋以改變Ω 角的自由度 (見 Fig.3-2)。我們可將這些轉動以一系列轉換矩陣來表示( , , )θ φ Ω →( ', ', ') θ φ Ω = ( +θ Δθ φ φ, +Δ Ω ΔΩ 變化： , + )

cos ' - sin ' 0 1 0 0 cos ' sin ' 0 cos( + ) -sin( + ) 0 1 0 0 = sin ' cos ' 0 0 cos ' sin ' - sin ' cos ' 0 sin( + ) cos( + ) 0 0 co

0 0 1 0 0 1 0 - sin ' cos ' 0 0 1 (3) (6) (5) (4) J φ φ φ φ φ φ φ φ θ θ φ φ φ φ θ θ ΔΩ ΔΩ ΔΩ ΔΩ ⎛ ⎞⎛ ⎞⎛ ⎞⎛ ⎞ ⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟⎜ _⎟⎝ _⎠ ⎝₁₄₄₂₄₄_{3 14}⎠⎝ ₄₂₄₄_{3 14}⎠⎝ ₄₂₄₄₃_{⎠ 144424443} 0 0 0 cos sin 0 s - sin -sin cos 0

0 sin cos 0 0 1 (2) (1) ' ' ' x x y J y z z φ φ θ θ φ φ θ θ ⎛ ⎞⎛ ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎝ ⎠⎝ ⎠ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ 1442443 14243，且。從Fig. 3-3 上之初始位置( , , )x₀ y₀ z₀ 開始，至Fig. 3-9 上之末位置( ', ', ')x y z 。我 們將以各示意圖說明每個轉換矩陣的功能：

(24)

Fig. 3-3 螺旋上一點( , , )x₀ y₀ z₀ 作轉動前之初使位置圖。其座標角度為( , , )θ φ Ω 。

x y z

為系統之座標軸，

x

h

y

h

z

h則為螺旋截面及自轉軸之座標軸。 Fig. 3-4 矩陣(1)作用示意圖。將螺旋之自轉軸

z

h，以順時針轉至ＹＺ平面上。同時點座標變為 1 1 1 ( , , )x y z ，且角度 (θ, 0, −Ω φ )。 Fig. 3-5 矩陣(2)作用示意圖。將螺旋自轉軸

z

h轉至平行Ｚ軸之方向。點座標變為( , , )x₂ y₂ z₂ ，且角度為( 0, 0, Ω−φ )。

z

x

y

1 1 1 ( , , )x y z φ Ω− θ h y h z h x

z

x

y

2 2 2 ( , , )x y z φ Ω− h y h x

(25)

Fig. 3-6 矩陣(3)作用示意圖。將螺旋繞著Ｚ軸逆時針轉(ΔΩ +φ) 度。使螺旋自轉面轉移到任意選定的角度( ' Ω 下，亦同時補償矩) 陣(1)所會影響到的自轉角變化 (Ω − + ΔΩ +φ φ = ')Ω 。點座標變為 3 3 3 ( , , )x y z ，且角度( 0, 0, Ω')。 Fig. 3-7 矩陣(4)作用示意圖。將螺旋繞著Ｚ軸順時針轉φ' 度，以預先補償稍後矩陣(6)所會影響之自轉角變化。點座標變為( , , )x₄ y₄ z₄ ，同時角度變為 ' ' ( 0, 0, Ω −φ )。 Fig. 3-8 矩陣(5)作用示意圖。將螺旋自轉軸

z

h 轉至任意選定的傾斜角( ' θ )下。點座標變為 5 5 5 ( , , )x y z ，同時角度變為 ' ' ( 'θ, 0, Ω −φ )。

z

x

y

3 3 3 ( , , )x y z h x ' Ω h y

z

x

y

4 4 4 ( , , )x y z h x ' φ' Ω − h y

x

y

5 5 5 ( , , )x y z ' φ' Ω− ' θ h y h z h x

z

(26)

Fig. 3-9 矩陣(6)作用示意圖。將螺旋之自轉軸

z

h，以逆時針轉至任意選定的方位角_{( )}φ ' _下。點座標變為轉動後之末位置 ( ', ', ')x y z ，角度成為( ' ' )θ φ, , Ω'。在被選定的螺旋移轉到新的位置後，系統便開始計算座標改變所造成的能量變化，再經由執行Metropolis 準則來判斷新結構是否將取代舊結構。因此在第二階段模擬中最重要的部份即是系統自由能的計算，我們必須將 RPs 在膜中摺疊的決定因素盡可能的剖析並將它一一量化，再來便是決定各個因素間的比重大

小。在我們先前的研究中 (Chen & Chen,2003;Chen et al.,2008)，已經成功地發展出幾項對蛋白質折疊模擬具重大影響性的能量項，再加上於本模型中新增之能量

項，都將於下面的段落中一一說明。

根據前人的研究指出 (Popot & Engelman,2000)，膜蛋白的穩定性可由兩個面向來考慮：一是螺旋與環境間的作用；二是螺旋間的作用。由於氫鍵、共價鍵或

離子對等較強的鍵結，並不常在螺旋間的作用裡被發現，所以影響螺旋間作用的

因素有兩個，一個是凡德瓦爾氏 (van der Waals，以下簡稱 vdW) 力，另一則是 R 基在不同構象 (rotamers) 間所造成的能量差異。因此我們所探討的第一個能量項便是不同螺旋上殘基間距離所造成的位勢能。vdW 位勢能可以利用 6-12 位能函數來描述： 1 0 0 0 0 12 6 , , (3-1)

r + r

E

=

{[

]

[

] }

r( , )

i j i j m n m n vdW i j m n i j i j i j

e

m n

≠

−

∑∑

，

x

y

' Ω ' θ h y h z

z

( ', ', ')x y z h x ' φ

(27)

其中e 為 vdW 作用項的強度，₁ m 與 i n 則分別代表第 i 和 j 根螺旋上的殘基序j 號。而在我們的模型中殘基位置即為相對應 Cα原子的位置，r

(

m ni, j

)

便代表上述兩殘基Cα間的距離 (如 Fig. 3-10 所示)。Tab. 3-1 中列出了 20 種基本胺基酸的 vdW 體積，其是以胺基酸內各原子的總數與其鍵結方式估算而出 (Nolting,2006:7;Richards,1974)。在模型中，我們將各殘基視為一具其對應 vdW 體積的球體，並以Cα原子作為圓心，故球體之半徑即為vdW 半徑。Eq. 3-1 中的 0 r i m 及 0 r j n 即代表其殘基所對應之vdW 半徑，而 0 0 (r + r ) i j m n 　便為兩殘基間，吸引力將等於斥力時的距離。若實際距離小於此值將會使其12 次項遠大於另一 6 次項，而使能量大幅升高。最後 Eq.3-1 為所有螺旋間的全部殘基組合計算加總而成，螺旋內部則因為 Cα 原子們相對位置保持不變而略去其間的 vdW 能量項不做計算。因為r 與0 r

(

m ni, j

)

皆將因螺旋鏈上胺基酸序列的不同而有所差異，序列相依的因子也在此被引進我們的模型中。 Fig. 3-10 vdW 位勢能作用示意圖。因為穿膜蛋白會接觸到的環境將有水的區域與雙層膜區域兩種 (如 Fig.3-2 所示)，我們將以兩個能量項來單獨討論螺旋與他們兩者之間的作用效應。首先是螺旋－水分子的作用。在先前許多以胜肽 (peptide) 與細胞膜系統所做的模型在模擬此效應時，便以一平滑介面的概念來區隔水與脂質兩種環境，再以厭水性

指數(hydropathy index) 來估算胜肽暴露在水中時的能量 (Biggin et al., 1997; La Rocca et al., 1999; Milik & Skolnick, 1993)，其中很常見的 Kyte-Doolittle 厭水性指數 (Tab.3-1)，是以各胺基酸在真空環境與水之間轉移時所需的 Gibbs 自由能，來量化其在水中時的厭水性 (hydrophobicity)，而其在脂質膜區域內的厭水性指 0 r i m　 i m　 ni 0 r j n r( , )m n_i _j

(28)

數則訂為零。我們採用了歸一化後的Kyte-Doolittle 指數 ( Tab.3-1 的 hy )來計算螺旋－水作用下的E_hw，其表示式為： w 2 , (3-2)

E =

hy

i h m i i m

e

∑

當殘基位於水中時

， 2 e 為此螺旋－水作用的強度，Eq. 3-2 便為將所有暴露在水中的殘基其對應的歸 一化後指數hy 加總再乘以一係數後的能量值。 Tab. 3-1 胺基酸性質表。上下順序以厭水性的由弱至強排列，厭水性指數採 Kyte-Doolittle 量表 (Kyte & Doolittle,1982)。

# _{hy 為模型所採用的指數量值，為 Kyte-Doolittle 歸一化後的結果，範圍由-1 到 1。} *資料來源 (Nolting,2006:7;Richards,1974)。vdW 半徑由 0 3 4 (r ) vdW = 3 π 體積推算。性質胺基酸 Kyte - Doolittle 厭水性指數 # hy (歸一化後指數) 3 * vdW 體積 (Αo ) 0_{( )} vdW r Αo 半徑 Arg ( R ) - 4.50 -1.00 148.00 3.28 Lys ( K ) - 3.90 -0.87 135.00 3.18 Asn ( N ) - 3.50 -0.78 96.00 2.84 Asp ( D ) - 3.50 -0.78 91.00 2.79 Gln ( Q ) - 3.50 -0.78 114.00 3.01 Glu ( E ) - 3.50 -0.78 109.00 2.96 His ( H ) - 3.20 -0.71 118.00 3.04 Pro ( P ) - 1.60 -0.36 90.00 2.78 Tyr (Y ) - 1.30 -0.29 141.00 3.23 Trp ( W ) - 0.90 -0.20 163.00 3.39 Ser ( S ) - 0.80 -0.18 73.00 2.59 Thr ( T ) - 0.70 -0.16 93.00 2.81 Gly ( G ) - 0.40 -0.09 48.00 2.25 Ala ( A ) 1.80 0.40 67.00 2.52 Met ( M ) 1.90 0.42 124.00 3.09 Cys ( C ) 2.50 0.56 86.00 2.74 Phe ( F ) 2.80 0.62 135.00 3.18 Leu ( L ) 3.80 0.84 124.00 3.09 Val ( V ) 4.20 0.93 105.00 2.93 Ile ( I ) 4.50 1.00 124.00 3.09

(29)

在接下來的螺旋－脂質膜的作用裡， Kessel 和 Ben-Tal 在綜合了眾多實驗與理論模擬的數據資料後指出 (Kessel & Ben-Tal,2002;Kessel et al.,2003)，當相對於脂質較為固態的穿膜螺旋存在於雙層膜中時，將使週遭的環境也變得較為“堅硬＂。即降低了系統的自由度 (或是熵)，而使 Gibbs 自由能隨之升高，增加的量將正比於脂質與螺旋間的接觸面積多寡。因此我們使用對應於θ 角的函數來量化此作用量，其能量函數表示為： 3 (3-3)

E =

h l

(1- cos )

i i

e

∑

θ

_，其中e 為此能量項的強度，₃ θ_i則為第 i 根螺旋的傾斜角，當θ_i越大時將使更多的螺旋表面區域埋浸在膜中，而讓E_{h l}連帶升高。接下來我們嘗試建構一綜合性描述螺旋－螺旋間和螺旋－脂質間作用的能

量。從文獻中得知 (McLean et al.,1991;Sugihara et al.,1996)，決定螺旋鏈在脂質

膜中構形的因素主要有三：(1)螺旋上非極性的面與脂質間的厭水性作用 (螺旋－脂質)；(2) 不同螺旋間極性端與極性端的親水性作用 (螺旋－螺旋)；(3)親水性殘基與水分子間的作用。其中第三項以包含於Eq.3-2 的厭水能量中。(1)與(2)則來自於螺旋在膜中可能同時出現厭水性與親水性作用的特質，一般將之稱為穿膜蛋白的兩性 (amphiphilic) α 螺旋結構特性，其作用與特定的胺基酸種類分布有關。至於胺基酸如何與為何如此分布，因近代蛋白質結構解析度的提升下，開始在文獻中被詳細地討論。其中提到了一些 R 基較大且厭水的胺基酸，例如 leucine、isoleucine、valine 和 phenylalanine，相對於埋在蛋白質的內部，他們較喜歡也較常與脂質接觸。其他具較小R 基的胺基酸，如 glycine、alanine、serine 和 threonine 則是頻繁地出現在膜蛋白中不同螺旋的介面裡 (Ulmschneider et al.,2005; Liang et al.,2005;Ben-Tal & Fleishman,2002)。這說明了胺基酸的種類對於穿膜蛋白的穩定結構具有一定的影響力。Ulmschneider 等人在他們的研究中 (2005)，同時也利用厭水性指數，來探討殘基在穿膜蛋白表面(或外部)與內部，

(30)

分布種類之不同所造成的厭水性表現差異。他們發現，即使存在於膜內的殘基其平均之厭水性指數，必定將比存在於水中之殘基之高，但就都處於膜內的殘基來說，尚發現他們將應映於蛋白質內不同的位置，而有親水性與厭水性之分布區別。且整體來說，分布在蛋白質表面的殘基其厭水性表現，將遠比分布在內部的殘基強得多。根據此穿膜蛋白兩性分布的特性與影響性，我們希望找到能將它放入能量計算中的物理量形式。首先分析待研究的兩個穿膜蛋白，BRD 與 HRD 的胺基酸序列對應之厭水性指數分布。我們發現在每根螺旋上，除了兩端容易探出膜外而暴露在水中之殘基以外，中間區段的殘基指數分布的確有一定的循環。配合我們利用ME 得到的各 α 螺旋結構，可發現厭水性指數為負值的殘基傾向於集中在螺旋結構的某面上。如果我們將所有螺旋的親水面 (具有較多親水性殘基分布的面)，都轉向 181°－ 270°的Ω 角象限，然後統計各螺旋上每個 Ω 角象限內 (1°－90°; 91°－180°; 181° －270°; 271°－360°)，厭水指數為負值之殘基出現的比例多寡，結果如 Fig.3-11 所示。將得到的比例結果與兩蛋白質 PDB 的原始結構圖，實際對照後發現，那些聚集較多親水性殘基的螺旋表面，都大致朝向蛋白質的內部，即通道的內裡部分。它的反面，也就是位於膜蛋白結構的表面位置，則有最高的厭水性殘基聚集密度。預測這樣特殊且一致的排列，將對了解穿膜蛋白內螺旋的Ω 角分布有極大 Fig. 3-11 親水性殘基出現比例與 α 螺旋上四個兩性特質面 (AFs) 投影示意圖。甜甜圈圖樣部分為Cα原子投影於ＸＹ平面時的情形，顆數便代表該面的殘基數。其上％數代表，BRD 與 HRD 兩蛋白質共 14 根 α 螺旋上，厭水性指數 < 0 的殘基在該象限內出現的比例。

(31)

的幫助。因此我們考慮依照螺旋表面的厭水或親水之特性，來建構不同螺旋表面

間的兩性交互作用 (amphiphilic interaction)。

藉著將殘基分成位於四個不同象限範圍內的想法，我們將每根α 螺旋的表面

分成4 個兩性特質面 (amphiphilic faces，以下簡稱 AFs)，如 Fig.3-11 所示。首先

扣除分布在螺旋兩端且具親水性的殘基之後，每個AF 上有最多六個，但平均為

五個的殘基數分布其中。我們再依每個 AF 上的殘基以對應的 hy 指數，沿Ｚ軸

位置由高至低，排列出每面的amphiphilic indexes。為了往後的兩性交互作用能

量計算之便，我們將擁有六個殘基數的組別，依照螺旋頭尾對此作用影響較弱的

原則下 (因較易離開雙層膜而進入水中，便不再是本能量項之討論範圍內)，採用了Eq. 3-4 將每面的 amphiphilic index 都降為五組：

1 2 2 3 3 4 4 5 5 6 0.5 hy +0.5 hy 0.5 hy +0.5 hy 0.5 hy +0.5 hy (3-4) 0.5 hy +0.5 hy 0.5 hy +0.5 hy × × ⎡ ⎤ ⎢ _× _× ⎥ ⎢ ⎥ ⎢ × × ⎥ ⎢ _× _× ⎥ ⎢ ⎥ ⎢ _× _× ⎥ ⎣ ⎦ 　　。試舉一例說明，以BRD 的 H2 上第 1 個 AF 而言，出現在上面的胺基酸順序由高至低是DKALATSL。扣掉了 hy 分別是 -0.78、-0.87 的 D 和 K 後，依照上式得： A L L A A T T S S L hy hy hy hy hy hy hy hy hy hy 0.5× + 0.5× 0.5× (+0.40) + 0.5× (+0.84) ⎡ ⎤ ⎢_0.5× _{+ 0.5×} ⎥ _{0.5× (+0.84) + 0.5× (+0.40)} ⎢ ⎥ ⎢0.5× + 0.5× ⎥ = 0.5× (+0.40) + 0.5× (−0.16) ⎢_0.5× _{+ 0.5×} ⎥ _{0.5× (−0.16) + 0.5× (−0.18)} ⎢ ⎥ ⎢_0.5× _{+ 0.5×} ⎥ _0.5 ⎣ ⎦ 　 +0.62 ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢_+0.62⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢+0.12⎥ ⎢ ⎥ ⎢_−0.17⎥ ⎢ ⎥ ⎢ ⎥ ⎢ _{× (−0.18) + 0.5× (+0.84)}⎥ ⎢_+0.33⎥ ⎣ ⎦ ⎣ ⎦ 。

此後，第 i 根螺旋上，第 j 個 AF 的第k位之amphiphilic index 表示式便為 amp_{, ,} i j k

ε

。各蛋白中的詳細分組情形我們將在本章第三節與第四節中說明。另外，各 AFs 的中心位置便是將其上所有殘基的Cα座標平均而得。由Ulmschneider 等人(2005) 的文獻中可得知，此兩性交互作用不僅是存在於螺旋之間，尚有螺旋－脂質的部份，且對於較為厭水的殘基來說，脂質的環境將更具有吸引力。因此相對於螺旋

(32)

上

ε

_{i j k}amp_{, ,} 的值將介於 -1 到 1 之間，我們將膜的 amphiphilic index 設為 1，表示式為

ε

_lipidamp。有了上述的物理量之後，接著我們建立不同螺旋間兩性特質面的交互作用。先以不同螺旋間如何可以成立接觸並發生作用的判斷條件開始，螺旋間的 AFs 相對位置如Fig. 3-12 所示。我們為此設下了兩個限制條件：(1) 兩螺旋的中心點 (各螺旋上 Cα原子座標之平均) 距離需滿足：r

(

i i₁, ₂

)

<30 A (3-5) o 。此一距離限制設得相當鬆散，以避免其過於干擾vdW 作用項的表現。況且 E_vdW 自然會排除螺旋間過遠或過近的距離。(2) 螺旋之間必須沒有被它根螺旋遮蔽的現象。在蛋白質摺疊的過程中或結果裡，同時有多根 α 螺旋彼此距離相近，如 Fig.3-12 中的 Hi - H₁ i (H₃ i 為第₁ i 根螺旋) 位置一般。我們以此二螺旋配對的中₁ 點夾角Δ ，來作為(2)之判斷條件。在考量了 α 螺旋的平均面徑並嘗試了幾次模ϕ 擬之後，我們發現當滿足： Δϕ ≥π 10 (3-6)， Hi 與 H₁ i 之間的接觸便沒有被 H₂ i 遮蔽的疑慮。 ₃ Fig. 3-12 α 螺旋截面圖與其上的 AFs 分布示意圖。三個空心大圓代表三根 α 螺旋： 1 2 3 H , H i i 與 Hi ，投影在ＸＹ平面上時的截面。四個實心小圓則為各螺旋上四組 AFs 的中心點。

(33)

當Eq. 3-5 和 Eq. 3-6 同時被滿足時，便可來篩選 Hi 與 H i₁ _２上將產生兩性交互作用之配對AFs。首先，Hi 上與 H₁ i 發生作用之 AF，將會是距離 H₂ i 中心點₂ 最近的 1 1 AF j i (Hi 上的兩性特質面1 j ，其上共有1 j -1 j )。因為兩 AFs 間的相互作4 用大小將與其相對方位有關，故我們將在 Fig. 3-12 中的 _AFj1 i₂ 至 AF 4 j i₂ 間進行角度篩選。雖說Hi 有₂ _AFj1 i₂ 至 AF 4 j i₂ 四種可能的作用面，但Fig. 3-12 中灰色的半圓區塊將因被遮蔽而無法與 1 1 AF j i 發生接觸，因此實際上應只有圓中白色的半圓區塊內的 AF，可以合理的和 1 1 AF j i 互相作用。故若 AF 11 j i 的自轉角為 = j 0 i11 Ω ，允許範圍內的 AF j i₂ 自轉角將會是 3 2 2 j i π ≥Ω ₂ ≥π 。換言之即是，若定義 1 , ij , ij = ij ij 1 1 1 2 2 ΔΩ Ω −Ω ，Hi 上可和₂ 1 1 AF j i 作用的面需符合： , ij , ij - π π 2 (3-7) 1 1 2 ΔΩ ≤ 　　　　　　　　。同時滿足此條件的有_AFj1 i₂ 和 AF j i₂2，我們的模型將選擇距離最近、角度也最吻合的_AFj1 i₂ ，來作為 Hi 上與 H2 i 發生作用的兩性特質面。這時我們稱1 AF ₁1 j i 與 1 AF j

i₂ 互為對方的可作用之兩性特質面 (interactive amphiphilic faces)。在我們的

模型中，此兩性交互作用項就如同 vdW 作用項一樣，每根螺旋將被容許與多根

螺旋同時發生反應。故就Fig.3-12 中的例子裡，因 Hi 也將滿足 Eq. 3-5 和 Eq. 3-6₃

的限制，其也構成與 Hi 作用之條件，模型便將依循挑選距離最近的和 Eq. 3-7₁ 中的原則來篩選其作用面，而很顯然地將會是 AFj i₃1與 AF ₁1 j i 發生作用。在所有螺旋都經過Eq.3-5 至 Eq.3-7 的條件判斷後，各螺旋上必然還是會有沒找到任何配對的面。其無法與其他螺旋產生接觸，代表此AF 必曝露在脂質的環境中，故將之歸類於螺旋－脂質的作用區塊。俱備了螺旋－螺旋與螺旋－脂質

(34)

(

)

1 , , 5 , , , , , , , , 1 4

+

E

=

(

)

amp i j k j j amp amp i i i i i j k i j k i i k i i amp

e

f

2

α

1 2

ε

1

ε

2

ε

1 2 1 2 = ≠

⎧⎡

⎤

⎛

⎞

⎪⎢

_ΔΩ

_Γ

⎥

⎨

_⎢

⎜

⎟

_⎥

⎝

⎠

⎪⎢

_⎣

⎥

_⎦

⎩

⋅

−

_∑

(

)

4 5 , , , 1 1 1 (3-8)

+

n

1 -

amp

i j i j k i j k

χ

ε

= = =

⎫

⎤

⎡

_⎪

⎥

⎢

⎬

⎥

⎢

⎪

⎣

∑ ∑

∑

_⎦⎭

，

其中e 為此能量項的強度。式中包含了兩個二維矩陣：₄ i, i n n F f_{1 2} × ⎡ ⎤ = ⎣ ⎦ 與 , i j ₄_n X

χ

× ⎡ ⎤ = ⎣ ⎦ ，兩矩陣的初始值皆設為零。其中i 為螺旋之足標；₁ j 為 AF₁ 之足標。n為目標蛋白質的穿膜 α 螺旋總數。 f_i_{1 2}_,_i 為矩陣 F 內，i 與₁ i 是否有₂ 產生兩性交互作用的表示值。若有作用其值為1，否則即為 0，其判斷準則即為上面提到的Eq.3-5 與 Eq.3-6。χ_{i j}_, 則特別針對某 AFj i ，是否有與其他螺旋上任何的 AFs 結成配對的表示值。不管實際配對的數目多少，只要有一組以上，值便為1，反之則維持初始值 0。因此就 Fig.3-12 裡的例子則為：χ_i₁_,_j₁= 1，χ_i₂_,_j₁= 1 和 = χ_i₃_,_j₁ 1。從上述兩個矩陣可看出 Eq.3-8 的第一項為螺旋－螺旋作用，第二項則為螺旋－脂質作用。另外一對應於Ω 角的函數：

(

,

)

, , , 2 = 1- ( - ) j j j j i i i i

α

π

1 2 1 2 ⎛ ⎞ ΔΩ _⎜ ΔΩ _⎟ ⎝ × ⎠ ，為可作用之兩性特質面間的夾角參數。因 - 2 j , _,j i1 i2 π π ΔΩ ≤ 為必然條件 (見 Eq.3-7)，故α 函數值將永遠介於 0 到 1 之間，且隨 , , - j j i1 i2 π ΔΩ 減小而升高。此項的用意為驅使模型找到方向最為吻合的接觸面來進行配對。找到正確的AFs 配對後，再以兩者的

ε

_{i j k}amp_{, ,} 來計算能量。Eq.3-8 中主要的部份： 5 _{, ,} _{, ,} 1 amp amp i j k i j k k

ε

₁

ε

₂ =

⋅

∑

，是將 1 AF j i 與 AF j i2 上各自的 , , amp i j k

ε

序列，依於膜中的對應高低位置兩兩相乘再加總 (見 Fig. 3-13 中之左圖)。因此若是對應的兩 amp_{, ,} i j k

ε

異號，即不是一個穩定的接

(35)

觸作用，將使E_amp整體升高，若同號則反之。

除了最相鄰的

ε

_{i j k}amp_{, ,} 產生作用之外，我們也嘗試了不同的作用範圍。也就是 Eq.3-8 中 ( amp_{, ,} )

i j k ε

Γ 修正項的作用。我們稱之為amphiphilic energy factor (以下簡寫為f_amp)。若 ( amp_{, ,} ) = 0 i j k ε Γ ，則f_amp =

(

0.0,1.0,0.0　

)

，其接觸作用便如Fig. 3-13 左圖一般。如果將作用面上兩相鄰的

ε

_{i j k}amp_{, ,} 也納入作用中，即Fig. 3-13 之右圖。但因距離較遠與相對方位角度較大，勢必使其具相對較弱的強度。故能量係數項將改為

(

)

famp = 0.2,1.0,0.2　　，而 , , ₁ ₂ ₁ ₂ 5 4 , , -1 , , , , , , 1 2 1 0.2 ( amp) i j k

amp amp amp amp

i j k i j k i j k i j k k k ε

ε

₊ = = = ⎛ ⎞ Γ _⎜_⎜ _⎟_⎟ ⎝

∑

⋅

+

∑

⋅

⎠。我們同時也嘗試了f_amp =

(

0.1,0.2,0.3,1.0,0.3,0.2,0.1　　　　

)

的能量項，其原理相同，但 , , ( amp) i j k

ε

Γ 項相當冗長故在此省略。處理完螺旋－螺旋區塊後，再將剩下沒配對的 AFs 與脂質做反應。道理完全一樣，但是更為簡單。因為我們將膜內的環境視為均勻且無限的，所以在螺旋－脂質的接觸作用裡，不需要考慮角度或距離的差異，僅需要將AFs 上所有 amp_{, ,} i j k

ε

乘上

ε

_lipidamp再行加總。但因後者恰為1，故在 Eq. 3-8

中被省略。需注意的是在Eq.3-8 的強度e 前有一負號，故不論是在螺旋－螺旋 ₄

Fig. 3-13 可作用之兩性特質面作用方式簡單示意圖。圓柱體為 α 螺旋外部構造假想圖，其上的五個實心圓點，代表五個

ε

_{i j k}amp_{, ,} 發生作用的大約排列位置。

(36)

或螺旋－脂質的部分，當同號的

ε

amp接觸而發生作用時，將使能量降低，反之則升高。

最後一能量項係有關於 retinal 分子對 RPs 摺疊結構的影響。在我們在近期

的研究裡 (Chen & Chen,2009) 曾表明， RPs 的七根螺旋很容易因為 vdW 接觸而形成六角形的摺疊，其中一根螺旋被其他六根包圍。雖然此一結構具有最低的

vdW 能量(Kokubo and Okamoto,2004)，但 retinal 分子將因 vdW 中的斥力，被從 RPs 的通道內排擠出結構外，而曝露在脂質的環境中。在脂質內的 retinal 分子不但提高了脂質分子的自由能，也無法與水分子形成氫鍵鍵結以及與螺旋鏈有較為穩定的非共價性接觸作用。因此我們嘗試在模型裡，將retinal 分子放入結構中來模擬此部分的接觸作用。首先，我們將retinal 分子的構型簡化成一棍狀的形狀，如圖Fig.3-13 所示。再將 C13-C15，C10-C12，C7-C9和C1-C6依序地分為組別I，II， III 和 IV，並按照順序地排列在棍狀結構上。為求精確，我們從 PDB 內的平均 retinal 座標數據中，取出C C₁₁ ₁₄(即 C11至 C14的距離)、C C8 11和C C5 8 ，來作為

鄰近組別間的距離 (標於 Fig. 3-13 中)。因 retinal 與 Lys 有共價鍵結構，組別 I 另外與Lys 以一定的鍵長連結而起，並隨 Lys 在膜內移動。retinal 本身則擁有θ 與

φ兩個自由度。為了讓此結構不至與其他螺旋碰撞而產生重疊，我們再依各組內

Fig. 3-14 Retinal 的原始化學結構與模型簡化結構示意圖。其上阿拉伯數字代表化

學結構中主鏈上之 C 數，羅馬數字則代表簡化後的原子團組別數，並分別標上鄰

(37)

鍵結種類與其數目來計算其V_vdW 與r_vdW (Abraham et al.,2003)，茲列於 Tab. 3-2。接著，我們考慮retinal 分子與螺旋的接觸能量 (retinal contact energy)：

( )

5 1 (3-9)

E

_contact ni _i i

e

ε

r

=

= −

∑

Δ

，其中e 為接觸作用強度。₅

ε

( )

Δ

r

i 是計算 retinal 與螺旋間是否接觸的函數，其中 i r Δ 為H i 中心點與 retinal 上四組別 (I - IV) 間的最短距離。若 6Å<Δ <11Å 則ri

( )

r_i =1

ε

Δ ，反之則為0。其中需注意的是，Eq.3-9 的主要目的為模擬 retinal－脂質與retinal－水間的能量差異。因此其概念為若 retinal 與越多的螺旋產生非共價性的接觸作用，則retinal 越可能位於 RPs 的通道內，而非暴露在脂質的環境中。接觸的螺旋本身不同並不會造成差異，故此項和其他能量項有所不同的是其序列不相依性。 Tab. 3-2 retinal 上四原子團之 vdW 體積 ( V_vdW) 與 vdW 半徑 ( r_vdW) 列表。其中 3 VvdW = 4 (rπ vdW) 3，計算方法參考自 (Abraham et al.,2003)。組別 I 組別 II 組別 III 組別 IV 3 V_vdW (Αo ) 72.24 56.47 72.24 139.47 rvdW ( )Α o 2.58 2.35 2.58 3.22

(38)

根據蛋白質摺疊的熱力學假說 (也稱 Anfinsen’s dogma)，蛋白質的原始結構將會是，其摺疊下具有自由能整體最小值狀態的構型 (Anfinsen,1973)。而我們的

任務，便是對應於上述的理論部分，將我們模型中所有針對 RPs 的結構特性發

展出的能量項，分別計算出來再行加總：

w (3-10)

E

_tot

=

E

_vdW

+

E

_h

+

E

_hl

+

E

_amp

+

E

_contact

。

再配合第二章所闡述之演算方法，模擬至整體能量最小的狀態。值得注意的是，

我們的模型是針對 RPs 的摺疊所發展而成。在本篇論文中我們選定了此蛋白質

家族中的兩名成員：BRD 與 HRD。兩者因胺基酸序列的不同，將會造成模擬中二級結構與部份指數的改變。我們將其間的細節與差異陳列在本章的下兩節中。

(39)

3-3 目標蛋白質 (1) － Bacteriorhodopsin ( BRD ) 在我們所做的模擬裡，因有首次嘗試的新能量項Eamp，其中牽涉了幾種不同的f_amp設定，也就是Γ(εi j kamp, , )項的改變。在期望先釐清此能量項作用細節的前提下，我們將BRD 中的四根 α 螺旋先行挑出，分別為 H3、H4、H5 和 H6。因它們位於整個BRD 蛋白質結構的轉角處，結構相對其他組合而言較為穩定，正巧適合較小規模也較單純化的蛋白質折疊模擬 (Baker at al.,2006)。希望可以藉此更清楚的釐清f_amp對模型的影響性。在之後的敘述裡，我們將此BRD 內 4 個螺旋

組成的子區塊 (4-Helix subdomain)，稱為 BRD4，來與 7 個螺旋的全區塊 (7-Helix domain)之 BRD7 做為區別。以下資料將直接以 BRD7 形式列出。有關 BRD4 的資訊，直接自H3、H4、H5 和 H6 的部份擷取即可。 Tab. 3-3 BRD7 的七根穿膜 α 螺旋上之胺基酸序列，與其二級結構對應 PDB 內原始二級結構之RMSD (僅 Cα原子)列表。胺基酸序列 RMSD(Å) H1 EWIWLALGTALMGLGTLYFLVKGM 0.861 H2 DPDAKKFYAITTLVPAIAFTMYLSMLL 1.172 H3 WARYADWLFTTPLLLLDLALLVD 1.206 H4 DQGTILALVGADGIMIGTGLVGALT 0.731 H5 YSYRFVWWAISTAAMLYILYVLFFGFTSKA 1.561 H6 RPEVASTFKVLRNVTVVLWSAYPVVWLIG 1.114 H7 PLNIETLLFMVLDVSAKVGFGLIL 1.469

穿膜蛋白質之結構預測

第一章 緒論

第二章 演算方法

−

(

)

E

>

E

<

{

}

(

)

(

)

⋅

(

)

第三章 模型

x y z

x

y

z

z

z

z

x

y

z

x

y

z

z

x

y

z

x

y

x

y

z

z

r + r

r + r

E

=

{[

]

[

] }

r( , )

r( , )

e

m n

m n

−

∑∑

x

y

z

(

)

(

)

E =

hy

e

∑

E =

(1- cos )

e

∑

θ

ε

ε

ε

(

)

(

第一章緒論

第二章演算方法

第三章模型

_ΔΩ

_Γ

_⎢

_⎥

_⎣

_⎦

_∑

_∑

_⎪

_⎦⎭