魏彥杰 張慧玲 黃慶生

(中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院高性能計(jì)算技術(shù)研究中心深圳市高性能數(shù)據(jù)挖掘重點(diǎn)實(shí)驗(yàn)室 深圳 518055)

蛋白質(zhì)折疊計(jì)算機(jī)模擬研究進(jìn)展

魏彥杰 張慧玲 黃慶生

(中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院高性能計(jì)算技術(shù)研究中心深圳市高性能數(shù)據(jù)挖掘重點(diǎn)實(shí)驗(yàn)室 深圳 518055)

蛋白質(zhì)折疊過程中的結(jié)構(gòu)變異可能導(dǎo)致“折疊病”，比如老年癡呆癥和多聚谷氨酰胺疾病等。因此蛋白質(zhì)折疊研究對(duì)于揭示“折疊病”致病機(jī)理、指導(dǎo)藥物設(shè)計(jì)等具有重大意義。文章闡述了蛋白質(zhì)折疊計(jì)算機(jī)模擬研究的研究近況，分別介紹了蛋白質(zhì)側(cè)鏈研究、蛋白質(zhì)折疊算法、蛋白質(zhì)折疊病研究、蛋白質(zhì)的分子動(dòng)力學(xué)模擬和蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)等幾個(gè)方面。

蛋白質(zhì)折疊；并行蒙特卡洛算法；蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)；分子動(dòng)力學(xué)模擬；蛋白質(zhì)折疊病

1 引 言

蛋白質(zhì)折疊問題是研究蛋白質(zhì)如何在短時(shí)間內(nèi)從一級(jí)結(jié)構(gòu)(一維多肽鏈)折疊為天然三維結(jié)構(gòu)，形成具有生命功能的大分子。生物體的遺傳信息(DNA)通過 RNA 轉(zhuǎn)錄和翻譯過程傳遞給蛋白質(zhì)，因此蛋白質(zhì)折疊也被稱為第二遺傳密碼，相關(guān)研究可以幫助揭示生命遺傳信息的表達(dá)和功能傳遞的奧秘[1]。在從一級(jí)結(jié)構(gòu)到天然三維結(jié)構(gòu)的折疊過程中，蛋白質(zhì)可發(fā)生誤折或聚集，其結(jié)構(gòu)和功能因此受到破壞，從而引起“折疊病”，比如老年癡呆癥、帕金森氏病和多聚谷氨酰胺疾病等[2]。美國(guó) 500 多萬(wàn)老年癡呆癥患者每年的花費(fèi)是 1830 億美元[3]，而中國(guó)的老年癡呆癥患者也有 500～600 萬(wàn)[4]。蛋白質(zhì)折疊研究對(duì)探索“折疊病”的致病機(jī)理意義重大，且有助于蛋白質(zhì)分子藥物的設(shè)計(jì)，故對(duì)“折疊病”的預(yù)防和治療也將起到重大的幫助作用[5]。因此我國(guó)將蛋白質(zhì)研究作為基礎(chǔ)科學(xué)的四個(gè)重大科學(xué)計(jì)劃之一列入“國(guó)家中長(zhǎng)期科學(xué)與技術(shù)發(fā)展規(guī)劃(2006-2020)”。

研究蛋白質(zhì)折疊的實(shí)驗(yàn)方法有 X 光晶體衍射，NMR 核磁共振等。雖然實(shí)驗(yàn)學(xué)方法積累了很多數(shù)據(jù)、貢獻(xiàn)巨大[6,7]，然而其耗時(shí)長(zhǎng)、費(fèi)用高，而且對(duì)于較難結(jié)晶的膜蛋白，很難得到其三維結(jié)構(gòu)[6-8]。另一方面，作為一個(gè)科學(xué)發(fā)現(xiàn)的獨(dú)立手段以及實(shí)驗(yàn)方法的重要補(bǔ)充，計(jì)算科學(xué)在蛋白質(zhì)折疊研究領(lǐng)域也起著越來(lái)越大的作用。近十幾年蛋白質(zhì)折疊的計(jì)算機(jī)研究方法取得了巨大發(fā)展，比如 Folding@home 利用世界范圍內(nèi)的計(jì)算機(jī)實(shí)現(xiàn)了蛋白質(zhì)折疊的分布式計(jì)算[9]，David Baker 研究組用蛋白質(zhì)折疊原理得到了全新構(gòu)象的蛋白質(zhì)[10]。本文將從側(cè)鏈對(duì)折疊的影響、高效的結(jié)構(gòu)空間搜索算法、折疊病致病機(jī)理、蛋白質(zhì)大體系的分子動(dòng)力學(xué)模擬以及蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)等領(lǐng)域介紹蛋白質(zhì)折疊的研究進(jìn)展。

2 蛋白質(zhì)側(cè)鏈研究

因?yàn)?20 種天然氨基酸的唯一區(qū)別是側(cè)鏈，所以側(cè)鏈對(duì)蛋白質(zhì)的結(jié)構(gòu)和功能，以及熱力學(xué)和動(dòng)力學(xué)特性非常重要[1]。蛋白質(zhì)折疊的側(cè)鏈研究主要是研究側(cè)鏈在蛋白質(zhì)折疊過程中所起的作用。

比較常用的蛋白質(zhì)側(cè)鏈的研究方法有以下幾種：(1)研究蛋白質(zhì)折疊中側(cè)鏈與側(cè)鏈之間氫鍵、雙硫鍵的形成[11]；(2)通過分割能量力場(chǎng)，研究側(cè)鏈與側(cè)鏈之間能量、側(cè)鏈與主鏈之間能量以及主鏈與主鏈之間的能量變化[12]；(3)研究側(cè)鏈熵變?cè)诘鞍踪|(zhì)折疊過程中的作用[13]；(4)利用已有的數(shù)據(jù)和貝葉斯統(tǒng)計(jì)方法等研究蛋白質(zhì)側(cè)鏈的統(tǒng)計(jì)規(guī)律，建立一個(gè)側(cè)鏈能量力場(chǎng)或者打分函數(shù)等[14]；(5)基于圓周統(tǒng)計(jì)方法研究側(cè)鏈的熱動(dòng)力學(xué)漲落現(xiàn)象[15]。目前這些研究方法在全面化和系統(tǒng)化方面尚有欠缺。其中，基于氫鍵或者雙硫鍵的方法要求氨基酸側(cè)鏈必需可以形成氫鍵或者雙硫鍵；基于能量分割的方法需要對(duì)能量力場(chǎng)的計(jì)算做出修改；基于熵的方法以及基于統(tǒng)計(jì)學(xué)(貝葉斯統(tǒng)計(jì)等)的方法需要對(duì)蛋白質(zhì)側(cè)鏈的結(jié)構(gòu)空間離散化處理，從而實(shí)現(xiàn)對(duì)不同狀態(tài)的側(cè)鏈結(jié)構(gòu)進(jìn)行統(tǒng)計(jì)。

蛋白質(zhì)計(jì)算機(jī)模型的構(gòu)建一般基于角度坐標(biāo)系或笛卡爾坐標(biāo)系。角度坐標(biāo)系中對(duì)側(cè)鏈的研究與笛卡爾坐標(biāo)系有所不同。由于角度具有2π 周期性，在角度坐標(biāo)系進(jìn)行蛋白質(zhì)折疊的計(jì)算機(jī)模擬時(shí)，對(duì)角度變量需要特殊處理。常用的方法是將角度或者角度差簡(jiǎn)單映射到 0～2π，或者 —π～+π[16]。但當(dāng)利用角度的平均值和方差研究蛋白質(zhì)的折疊過程時(shí)，基于線性統(tǒng)計(jì)學(xué)的方法會(huì)有許多問題，比如會(huì)導(dǎo)致一個(gè)單峰分布的角度函數(shù)變?yōu)殡p峰分布函數(shù)。如圖 1 所示，黑點(diǎn)和綠點(diǎn)在空間上本屬同一個(gè)分布，但由于角度范圍的選擇，黑點(diǎn)和綠點(diǎn)成為兩個(gè)不同的分布(兩個(gè)分布的角度平均值約為—160°和+160°)。

圖 1 線性統(tǒng)計(jì)導(dǎo)致錯(cuò)誤的數(shù)據(jù)分析Fig. 1. Analysis of circular data using linear statistics

相比于線性統(tǒng)計(jì)學(xué)，圓周統(tǒng)計(jì)學(xué)更適用于此類問題。但是圓周統(tǒng)計(jì)方差和協(xié)方差的數(shù)學(xué)模型在蛋白質(zhì)折疊研究中還存在不一致性[15]。圓周統(tǒng)計(jì)中角度 a 可以用二維向量表示為角度方差可以用來(lái)研究蛋白質(zhì)折疊過程中物理量的熱動(dòng)力學(xué)漲落，其公式為：

由此導(dǎo)出的方差為：兩個(gè)方差公式(1)和(3)不一致，因此圓周統(tǒng)計(jì)方差和協(xié)方差之間在研究蛋白質(zhì)的兩面角時(shí)存在不一致性，兩方差之間不存在線性統(tǒng)計(jì)中的對(duì)應(yīng)關(guān)系[15]。文獻(xiàn)[15]中研究角度波動(dòng)在較小的范圍內(nèi)，作者采用了公式(1)計(jì)算方差，用如下公式(4)來(lái)計(jì)算角度協(xié)方差：

上述計(jì)算方差和協(xié)方差的方法并不一致，且只在角度波動(dòng)較小的情況下適用。而波動(dòng)較大的角度對(duì)蛋白質(zhì)折疊的影響更大，因而具有更大的研究?jī)r(jià)值。

3 蛋白質(zhì)折疊空間搜索算法

對(duì)于一個(gè)只有 100 個(gè)氨基酸的蛋白質(zhì)，其構(gòu)象數(shù)目大約有 10100(假設(shè)每個(gè)氨基酸只有 10 個(gè)構(gòu)象)。 假設(shè)搜索一個(gè)構(gòu)象需要 10—15秒(使用目前最快的千萬(wàn)億次計(jì)算機(jī))，則搜索所有的蛋白質(zhì)結(jié)構(gòu)空間需約 3×1077年。因此搜索所有結(jié)構(gòu)空間不現(xiàn)實(shí)，這就對(duì)蛋白質(zhì)結(jié)構(gòu)高效搜索算法提出了更高的要求[17]。蛋白質(zhì)折疊算法依賴于一個(gè)基本假設(shè)，即蛋白質(zhì)的天然結(jié)構(gòu)是蛋白質(zhì)自由能最低的結(jié)構(gòu)[18]。而傳統(tǒng)的分子動(dòng)力學(xué)模擬和蒙特卡洛模擬方法在研究蛋白質(zhì)折疊時(shí)，低溫下的模擬都會(huì)被“陷”在蛋白質(zhì)能量的局部最優(yōu)區(qū)間，很難“跳”出來(lái)從而找到全局最優(yōu)解[19]。

廣義系綜方法(Generalized Ensemble)是解決此問題最常用的一類方法[19]。此方法的核心是使用非波爾茲曼分布函數(shù)實(shí)現(xiàn)模擬仿真在能量空間的自由行走，從而搜索更為廣泛的結(jié)構(gòu)空間，同時(shí)還可計(jì)算任一溫度下的正則系綜物理量。這類方法包括多正則系綜(Multicanonical)蒙特卡洛方法[20]、Wang-Landau 蒙特卡洛方法[21]、模擬回火(Simulated Tempering)[22]、1/k 搜索[23]和副本交換方法(Replica Exchange Method)[24]等。由于廣義系綜方法的非波爾茲曼分布函數(shù)未知，因此需設(shè)計(jì)一個(gè)繁瑣的程序決定此分布函數(shù)[19]。并行回火和 Wang-Landau 方法可較為容易地獲取非波爾茲曼分布函數(shù)。副本交換方法(又稱并行回火方法，Parallel Tempering)要求在多個(gè)溫度下分別進(jìn)行蒙特卡洛或者分子動(dòng)力學(xué)仿真，經(jīng)過一定時(shí)間后，不同溫度下的仿真以一定概率互換結(jié)構(gòu)信息；非波爾茲曼分布函數(shù)就是各個(gè)溫度下的波爾茲曼分布函數(shù)的乘積[19]。隨著蛋白質(zhì)系統(tǒng)的增大，此方法需要更多不同溫度下的仿真。而選擇多少仿真溫度，以及如何選取溫度成為此算法的一個(gè)難點(diǎn)。通過迭代修改一個(gè)算法參數(shù) F，Wang-Landau 算法可以自動(dòng)獲得非波爾茲曼分布函數(shù)以及系統(tǒng)狀態(tài)函數(shù)，這極大地簡(jiǎn)化了獲取非波爾茲曼分布函數(shù)的難度。然而此算法要求預(yù)先知道模擬系統(tǒng)的能量范圍，這樣才能實(shí)現(xiàn)模擬仿真在此能量范圍內(nèi)的自由行走[1]。

以上描述的算法大都在分子動(dòng)力學(xué)模擬和蒙特卡洛模擬中實(shí)現(xiàn)[19]。比如并行回火蒙特卡洛方法和并行回火分子動(dòng)力學(xué)模擬[25,26]。分子動(dòng)力學(xué)模擬可以研究系統(tǒng)的動(dòng)力學(xué)過程；而蒙特卡洛方法可以研究系統(tǒng)的熱動(dòng)力學(xué)過程。高精確度的蛋白質(zhì)折疊全原子模型模擬需要計(jì)算幾萬(wàn)到幾十萬(wàn)個(gè)原子之間的相互作用力，所以分子動(dòng)力學(xué)模擬一般只能模擬納秒級(jí)的折疊過程，因而對(duì)于研究發(fā)生在微妙到毫秒時(shí)間內(nèi)的蛋白質(zhì)折疊問題存在局限性。一般分子動(dòng)力學(xué)模擬從一個(gè)實(shí)驗(yàn)結(jié)構(gòu)開始模擬蛋白質(zhì)折疊過程[27]，而蒙特卡洛方法則可以研究蛋白質(zhì)折疊的整個(gè)熱動(dòng)力學(xué)過程，而且蒙特卡洛方法不依賴于初始結(jié)構(gòu)的選擇，可以搜索更廣泛的結(jié)構(gòu)空間。

其他用于蛋白質(zhì)折疊的研究算法有：結(jié)構(gòu)空間退火算法(CSA)[28]、遺傳算法[29]、 aBB[30]和圖論算法[31]等。其中，結(jié)構(gòu)空間退火算法首先搜索大量的可能含最優(yōu)解的子空間，然后在這些子空間中尋找最優(yōu)解[28]；遺傳算法的操作對(duì)象是一群蛋白質(zhì)結(jié)構(gòu)，通過蛋白質(zhì)結(jié)構(gòu)的突變、結(jié)構(gòu)域選擇和重組，實(shí)現(xiàn)蛋白質(zhì)結(jié)構(gòu)的進(jìn)化，最終找到蛋白質(zhì)的最優(yōu)結(jié)構(gòu)[28]；aBB 算法在理論上可以確保全局最優(yōu)解[29]。然而這些算法僅限于尋找蛋白質(zhì)的最優(yōu)解，不能研究蛋白質(zhì)折疊的動(dòng)力學(xué)和熱動(dòng)力學(xué)過程。

越來(lái)越多的蛋白質(zhì)折疊研究采用并行算法，比如并行回火算法[25,26]、CSA 和 aBB 混合并行算法[33]、并行遺傳算法等[33]。斯坦福大學(xué)Folding@home 項(xiàng)目利用世界范圍內(nèi)的計(jì)算機(jī)實(shí)現(xiàn)了蛋白質(zhì)折疊的分布式計(jì)算研究[9]。通過在世界各地的電腦運(yùn)行一個(gè)客戶端程序，F(xiàn)olding@ home 將世界上各地的電腦連在一起，構(gòu)成了世界上最大的超級(jí)計(jì)算機(jī)之一。

蒙特卡洛方法可以研究蛋白質(zhì)折疊的熱動(dòng)力學(xué)性質(zhì)，而且比分子動(dòng)力學(xué)模擬能搜索更廣泛的結(jié)構(gòu)空間。相反，分子動(dòng)力學(xué)模擬更適合研究蛋白質(zhì)的精細(xì)結(jié)構(gòu)變化，比如蛋白質(zhì)與化學(xué)小分子的相互作用。

另外一種研究蛋白質(zhì)折疊的思路是設(shè)計(jì)特殊的硬件設(shè)備，用于快速計(jì)算和模擬蛋白質(zhì)大分子的熱力學(xué)和動(dòng)力學(xué)特性。較早的例子是日本的 RIKEN研究機(jī)構(gòu)于 2006年開發(fā)的MDGRAPE-3，它主要用于特殊蛋白質(zhì)結(jié)構(gòu)的分子動(dòng)力學(xué)模擬。MDGRAPE-3 由 4824 個(gè)特制的 MDGRAPE-3 芯片以及英特爾 Xeon 處理器構(gòu)成，對(duì)分子的計(jì)算通過 MDGRAPE 芯片和通用芯片共同完成。在2006 年時(shí)，它已經(jīng)具備了千萬(wàn)億浮點(diǎn)運(yùn)算能力，比當(dāng)時(shí)超算 TOP 500 排名第一的 IBM Blue Gene/L 還要快 3 倍。對(duì)兩個(gè)不同的處于周期邊界條件下的生物分子系統(tǒng)的模擬表明，單個(gè) MDGRAPE-3 系統(tǒng)比傳統(tǒng)基于 CPU的系統(tǒng)要快 30～40 倍[34]。2009 年 Shaw 發(fā)表了面向分子動(dòng)力學(xué)模擬的專用計(jì)算機(jī) ANTON[35]。ANTON 包含 512 個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)是一個(gè)專用計(jì)算單元(ASIC)，包含一個(gè)高通量相互作用子系統(tǒng)和一個(gè)可變子系統(tǒng)，前者可以快速計(jì)算原子間的兩兩相互作用，后者則是可編程的子系統(tǒng)，用于完成負(fù)荷較輕的其他一些計(jì)算。512 個(gè)節(jié)點(diǎn)之間的總線互聯(lián)也考慮了分子動(dòng)力學(xué)模擬中空間區(qū)域分割的需求。ANTON 特殊的硬件配置給人留下深刻印象，但其實(shí)它是軟硬件緊密結(jié)合的一個(gè)平臺(tái)，設(shè)計(jì)者重新處理了分子動(dòng)力學(xué)模擬并行化時(shí)碰到的負(fù)載均衡問題，用了新的辦法進(jìn)行模擬的空間區(qū)域分割，甚至專門設(shè)計(jì)了定點(diǎn)算術(shù)運(yùn)算法則。與 MDGRAPE-3 不同的是，所有的ANTON 計(jì)算均在 ASICs 上完成。對(duì)于 23,558 原子構(gòu)成的分子系統(tǒng)的分子動(dòng)力學(xué)模擬，ANTON在一天內(nèi)可以模擬 17 微秒，而基于通用芯片的系統(tǒng)只能模擬幾百納秒。與 ANTON 相比，ANTON2 在軟件和硬件方面做了許多提升，其分子動(dòng)力學(xué)模擬的效率也因此提升了 16%[36]。

4 分子動(dòng)力學(xué)模擬

分子動(dòng)力學(xué)(Molecular Dynamics)模擬是一種數(shù)值模擬方法，通過將分子抽象為由化學(xué)鍵連接的質(zhì)點(diǎn)，按照基于牛頓力學(xué)的數(shù)學(xué)模型(力場(chǎng))迭代求解分子體系的行為。這種方法比量子力學(xué)方法大大節(jié)省了計(jì)算資源，并且可以求出分子體系隨時(shí)間演化的行為。在 20 世紀(jì)五、六十年代分子動(dòng)力學(xué)模擬被提出來(lái)時(shí)，它只能應(yīng)用于簡(jiǎn)單的單原子分子的同質(zhì)化體系，解決理論物理學(xué)的問題[37,38]。隨著它在化學(xué)和生物學(xué)中的廣泛應(yīng)用，模擬體系的原子數(shù)目越來(lái)越多、模擬時(shí)間越來(lái)越長(zhǎng)成為自然的需求，而計(jì)算機(jī)硬件和軟件技術(shù)的發(fā)展以及分子動(dòng)力學(xué)模擬算法本身的優(yōu)化則提供了這樣的可能性。在當(dāng)代，應(yīng)用于生物學(xué)的模擬體系原子數(shù)目通常在 10 萬(wàn)以上，可以是同時(shí)包含生物大分子、生物膜、溶劑和無(wú)機(jī)離子的復(fù)雜體系。

巨大的模擬體系要求投入大量的計(jì)算資源，為了在可接受的時(shí)間內(nèi)完成模擬，并行計(jì)算是現(xiàn)實(shí)的選擇。NAMD2 是并行化而且可伸縮性很強(qiáng)的分子動(dòng)力學(xué)模擬軟件。2013 年 Zhao等[39]在 Blue Waters 超級(jí)計(jì)算機(jī)的支持下，利用 NAMD2 軟件進(jìn)行了包含多達(dá)六千萬(wàn)個(gè)原子的模擬，協(xié)助冷凍電鏡解析出了人免疫缺陷病毒-1(Human Immunode fi ciency Virus-1，HIV-1)衣殼的高分辨率結(jié)構(gòu)。HIV-1 是一種逆轉(zhuǎn)錄病毒，成熟的病毒表面包被衣殼。衣殼由約 1300 個(gè)衣殼蛋白組成。衣殼蛋白相互連接，組成六聚物單元，六聚物單元又相互連接形成六聚物的五聚物(Pentamer-of-Hexamers，POH)或六聚物的六聚物(Hexamer-of-Hexamers，HOH)網(wǎng)格。HOH 網(wǎng)格是規(guī)整的長(zhǎng)管狀，衣殼大體上由 HOH 網(wǎng)格組成，而在特定的 12 個(gè)頂點(diǎn)嵌入 POH 網(wǎng)格，在這些頂點(diǎn)處衣殼的表面彎曲閉合。因?yàn)橐職さ恼w形狀不規(guī)則，缺乏對(duì)稱性，用冷凍電鏡三維重構(gòu)技術(shù)只能獲得分辨率為 8.6? 的結(jié)構(gòu)。Zhao 等采用了分子動(dòng)力學(xué)柔性擬合(Molecular Dynamics Flexible Fitting，MDFF)技術(shù)，將電鏡的電子密度圖作為外部勢(shì)，驅(qū)動(dòng)柔性的衣殼蛋白填充電子密度圖，得到 HOH 的原子分辨率結(jié)構(gòu)。根據(jù)由不同數(shù)目 HOH 和 POH 組成的幾種網(wǎng)格的形狀與衣殼的大體結(jié)構(gòu)比較，初步推出衣殼由 216個(gè) HOH 和 12 個(gè) POH 組成。為了得到完整的衣殼結(jié)構(gòu)，將 MDFF 得到的 HOH 結(jié)構(gòu)和之前報(bào)道的 POH 結(jié)構(gòu)組裝成成兩個(gè)模型并進(jìn)行分子動(dòng)力學(xué)模擬，其中由 216 個(gè) HOH 和 12 個(gè) POH 組成的衣殼的模擬體系共包含 64,332,531 個(gè)原子，由186 個(gè) HOH 和 12 個(gè) POH 組成的衣殼的模擬體系共包含 64,423,983 個(gè)原子，這些模擬各自持續(xù)了 100 ns。通過比較和分析，確認(rèn)了由 216 個(gè)HOH 和 12 個(gè) POH 組成的衣殼的原子分辨率結(jié)構(gòu)。成功模擬如此龐大的生物大分子體系標(biāo)志著分子動(dòng)力學(xué)模擬軟件和并行計(jì)算技術(shù)的成功結(jié)合。

雖然并行計(jì)算可以加快大體系的模擬速度，但分子動(dòng)力學(xué)模擬最根本的速度瓶頸在于通用CPU 串行地處理原子間的相互作用。要提升單機(jī)(單節(jié)點(diǎn))的計(jì)算速度，可以設(shè)計(jì)專用計(jì)算單元，挖掘計(jì)算機(jī)硬件的潛力。前文提到的面向分子動(dòng)力學(xué)模擬的專用計(jì)算機(jī) ANTON 就是最好的例子。ANTON 可以輕松地進(jìn)行微秒級(jí)的分子動(dòng)力學(xué)模擬，速度超越了同時(shí)代機(jī)器的兩個(gè)數(shù)量級(jí)。借助這個(gè)機(jī)器，Shaw 等[40]用長(zhǎng)達(dá) 1 ms 的分子動(dòng)力學(xué)模擬探討了蛋白質(zhì)折疊的問題；Arkhipov 進(jìn)行了長(zhǎng)達(dá) 200 μs 的模擬揭示了 EGF 受體和細(xì)胞膜的相互作用[41]；Dror 等[42]模擬了小分子配體和 G 蛋白偶聯(lián)受體自發(fā)結(jié)合的過程，研究了基于別構(gòu)效應(yīng)的藥物設(shè)計(jì)。

除了設(shè)計(jì)專門的軟件和硬件提升模擬速度，改進(jìn)力場(chǎng)和模擬算法本身也是重要的研究課題。由于力場(chǎng)和模擬算法的改進(jìn)，計(jì)算的速度、精確度和適用的分子體系都有很大進(jìn)步。常用于模擬生物學(xué)體系的 AMBER 力場(chǎng)和 CHARMM 力場(chǎng)都有兩類能量項(xiàng)：相鄰原子的鍵長(zhǎng)、鍵角和二面角，和原子兩兩相互作用的靜電力和范德華力。靜電力和范德華力是遠(yuǎn)程相互作用，因此要對(duì)體系中所有原子的兩兩組合求和，這對(duì)大體系是非常耗時(shí)的步驟。然而，遠(yuǎn)程相互作用的強(qiáng)度隨距離減弱，并且體系中的溶劑有一定的均一性，因此可以利用這些特點(diǎn)設(shè)計(jì)一些快速算法。Price和 Brooks 提出改進(jìn)的水溶劑模型 TIP3P，可以配合 Ewald 求和簡(jiǎn)化水分子中的氫原子的處理[43]。Cerutti 和 Case 設(shè)計(jì)了多水平 Ewald 求和算法，進(jìn)一步優(yōu)化了基于快速傅里葉變換的傳統(tǒng)算法[44]。另一方面，一直以來(lái)分子動(dòng)力學(xué)模擬方法無(wú)法處理共價(jià)鍵的斷裂和形成事件，主要是因?yàn)檫@些事件涉及化學(xué)變化，因此發(fā)生電子轉(zhuǎn)移，不符合分子動(dòng)力學(xué)模擬方法的基本假設(shè)。為了模擬更復(fù)雜的現(xiàn)象，Walker、Crowley 和 Case 給 Amber 補(bǔ)充了 QM/MM 的功能，兼顧了分子動(dòng)力學(xué)模擬方法的高效率和量子力學(xué)方法的靈活性[45]。

總之，分子動(dòng)力學(xué)模擬已經(jīng)是化學(xué)和生物學(xué)研究的實(shí)用技術(shù)，而這項(xiàng)技術(shù)也由于應(yīng)用的驅(qū)動(dòng)快速地發(fā)展，能越來(lái)越快地處理越來(lái)越大和越來(lái)越復(fù)雜的問題。

5 蛋白質(zhì)折疊與疾病

蛋白質(zhì)結(jié)構(gòu)決定了蛋白質(zhì)的功能，蛋白質(zhì)功能又與生理密切相關(guān)。蛋白質(zhì)折疊過程中的結(jié)構(gòu)變異可能導(dǎo)致“折疊病”，比如老年癡呆癥和多聚谷氨酰胺疾病等。因此對(duì)蛋白質(zhì)折疊和結(jié)構(gòu)的研究有助于研究許多疾病的致病機(jī)理，幫助藥物設(shè)計(jì)，從而找到治療方法。下面以多聚谷氨酰胺(polyQ)疾病為例介紹蛋白質(zhì)折疊與相關(guān)疾病的關(guān)系研究。

多聚谷氨酰胺疾病是因三核苷酸 CAG 基因異常重復(fù)擴(kuò)增導(dǎo)致產(chǎn)生多聚谷氨酰胺蛋白的神經(jīng)系統(tǒng)功能障礙疾病[46]。九種多聚谷氨酰胺疾病的共同特征是患者的大腦細(xì)胞中形成了不可溶的蛋白質(zhì)聚集體。雖然這九種疾病的致病蛋白質(zhì)以及相應(yīng)的基因各不相同，但是他們都有共同的多聚谷氨酰胺片段。至今尚未清楚多聚谷氨酰胺蛋白質(zhì)是如何誤疊和聚集，進(jìn)而導(dǎo)致疾病[47]?？赡艿闹虏C(jī)理如圖 2 表示，與其他許多神經(jīng)退行性疾病(比如老年癡呆證)相似，多聚谷氨酰胺疾病與多聚谷氨酰胺蛋白的誤疊有關(guān)[47]。當(dāng)多聚谷氨酰胺蛋白長(zhǎng)度達(dá)到 35～40 個(gè)氨基酸時(shí)，就容易產(chǎn)生這種疾病，而且多聚谷氨酰胺蛋白長(zhǎng)度越長(zhǎng)疾病就越嚴(yán)重[48]。雖然體內(nèi)和體外實(shí)驗(yàn)也支持蛋白質(zhì)的誤疊導(dǎo)致了疾病的假設(shè)，但各個(gè)實(shí)驗(yàn)的結(jié)果仍存在不一致的地方[49]。由于不可溶性，科學(xué)家還沒得到蛋白質(zhì)聚集的實(shí)驗(yàn)結(jié)構(gòu)[50]。

圖 2 多聚谷氨酰胺疾病可能致病機(jī)理Fig. 2. Possible mechanism of polyQ diseases

用非全原子分子動(dòng)力學(xué)模擬方法[48,51]發(fā)現(xiàn)當(dāng)多聚谷氨酰胺蛋白質(zhì)長(zhǎng)度超過 37 個(gè)氨基酸時(shí)，多聚谷氨酰胺蛋白質(zhì)形成 β 螺旋結(jié)構(gòu)；同時(shí)Khare1 等[48]指出側(cè)鏈和主鏈之間的氫鍵作用對(duì)于 β 螺旋結(jié)構(gòu)的形成起著非常重要的作用，而且β 螺旋結(jié)構(gòu)內(nèi)部結(jié)構(gòu)緊密，這區(qū)別于著名的內(nèi)部空間疏松的 β 螺旋模型[52]；Chopra 等[53]使用蒙特卡洛方法研究發(fā)現(xiàn)幾個(gè)關(guān)鍵的氨基酸決定了多聚谷氨酰胺蛋白質(zhì) β 螺旋結(jié)構(gòu)的形成，而且一個(gè)穩(wěn)定的 β 螺旋是蛋白質(zhì)聚集的基礎(chǔ)。

利用并行回火分子動(dòng)力學(xué)方法和非全原子模型，Nelson 等[54]研究了一個(gè)被廣泛應(yīng)用的蛋白質(zhì)系統(tǒng) CI2-polyQ，發(fā)現(xiàn)側(cè)鏈之間的氫鍵是蛋白質(zhì)聚集形成的關(guān)鍵。基于 cross-beta-spine 的空間拉鏈模體[55,56]建立了長(zhǎng)短不同的多聚谷氨酰胺蛋白模型；基于這些模型的分子動(dòng)力學(xué)模擬，Esposito[56]發(fā)現(xiàn)短的多聚谷氨酰胺也可以形成蛋白質(zhì)聚集，而且氫鍵對(duì)蛋白質(zhì)聚集的穩(wěn)定性作用巨大。

6 蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)

目前解析蛋白質(zhì)三維結(jié)構(gòu)的生物學(xué)實(shí)驗(yàn)方法主要有 X-RAY 和 NMR 法，但這些方法不僅復(fù)雜耗時(shí)，而且花費(fèi)較高。由于實(shí)驗(yàn)解析法存在這些不足，使得計(jì)算方法的發(fā)展成為必然。蛋白質(zhì)三維結(jié)構(gòu)預(yù)測(cè)的計(jì)算方法主要分為三類：同源模建法(Comparative Modeling)[57]、折疊識(shí)別法(Threading and Fold Recognition)[58,59]和從頭預(yù)測(cè)法(ab Initio Modeling)[60,61]。同源模建法根據(jù)相關(guān)蛋白質(zhì)的序列是否具有較高的相似性，判斷它們是否具有相似的三維結(jié)構(gòu)，預(yù)測(cè)的蛋白質(zhì)結(jié)構(gòu)具有很高的可信度。折疊識(shí)別法將未知蛋白的序列放入結(jié)構(gòu)數(shù)據(jù)庫(kù)中，并選擇最合適的折疊，以此來(lái)預(yù)測(cè)蛋白質(zhì)的三維結(jié)構(gòu)，但如果一個(gè)預(yù)測(cè)所需要的蛋白質(zhì)折疊不存在于折疊文庫(kù)中，該方法將失效。從頭預(yù)測(cè)法預(yù)測(cè)程序使用能量最小化原理，搜索每一種可能構(gòu)象來(lái)確定一個(gè)具有最低全局能量的構(gòu)象。

蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)技術(shù)評(píng)估(C r i t i c a l Assessment of Structure Prediction，CASP)大賽是一個(gè)世界性的蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)技術(shù)評(píng)比活動(dòng)，被譽(yù)為蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)領(lǐng)域的奧林匹克競(jìng)賽，自 1994 年開始，每?jī)赡杲M織一屆。盡管近年來(lái)在 CASP 中，不基于模板的預(yù)測(cè)方法(Free Modeling)在預(yù)測(cè)精度方面的提升不大，但由于被解析的蛋白質(zhì)結(jié)構(gòu)不斷增多，基于模板的預(yù)測(cè)方法(Template Based Modeling)在預(yù)測(cè)的準(zhǔn)確度上取得了較大的進(jìn)步[62]。大賽還將蛋白質(zhì)殘基作用關(guān)系預(yù)測(cè)(Contact Prediction)作為蛋白質(zhì)三維結(jié)構(gòu)預(yù)測(cè)比賽的一部分[63]。與往屆 CASP 相比，CASP 10 的“Re fi nement Category”中，Mirjalili等[64]使用了分子動(dòng)力學(xué)的方法對(duì)預(yù)測(cè)的蛋白質(zhì)結(jié)構(gòu)進(jìn)行優(yōu)化，并首次實(shí)現(xiàn)了對(duì)所有 CASP 10 目標(biāo)蛋白預(yù)測(cè)精度的整體提升。

7 總 結(jié)

本文從蛋白質(zhì)折疊中的側(cè)鏈研究、蛋白質(zhì)折疊算法研究、蛋白質(zhì)分子動(dòng)力學(xué)模擬以及蛋白質(zhì)折疊病研究等方面詳細(xì)介紹了蛋白質(zhì)折疊計(jì)算機(jī)模擬研究的進(jìn)展，這些進(jìn)展推動(dòng)了蛋白質(zhì)功能的研究和基于蛋白質(zhì)折疊的藥物設(shè)計(jì)，表明計(jì)算機(jī)模擬方法在生命科學(xué)的研究中扮演著越來(lái)越重要的角色。

[1] McKee T, McKee JR. 生物化學(xué)導(dǎo)論 [M]. 北京：科學(xué)出版社, McGraw-Hill出版社, 2000.

[2] Selkoe DJ. Folding proteins in fatal ways [J]. Nature, 2003, 426: 900-904.

[3] The Alzheimer’s Association. Alzheimer’s Disease Facts and Figures [EB/OL]. http://www.alz.org/ alzheimers_disease_facts_and_ fi gures.asp.

[4] 范東輝. 嚴(yán)重影響健康的老年癡呆癥 [J].中國(guó)檢驗(yàn)檢疫, 2006, 12: 64.

[5] Motta A, Reches M, Pappalardo L, et al. The preferred conformation of the tripeptide Ala-Phe-Ala in water is an inverse Gamma-turn: implications for protein folding and drug design [J]. Biochemistry, 2005, 44 (43): 14170-14178.

[6] Wuthrich K. Protein structure determination in solution by NMR spectroscopy [J]. The Journal of Biological Chemistry, 1990, 265(36): 22059-22062. [7] Drenth J. Principles of Protein X-Ray Crystallography [M]. New York: Springer-Verlag Inc., 1990.

[8] White SH. The progress of membrane protein structure determination [J]. Protein Science, 2004, 13: 1948-1949.

[9] Larson SM, Snow CD, Shirts M. Folding@ Home and Genome@ Home: Using distributed computing to tackle previously intractable problems in computational biology [J]. Proceedings of the Computational Genomics, 2002.

[10] Kuhlman B, Dantas G, Ireton GC, et al. Design of a novel globular protein fold with atomic-level accuracy [J]. Science, 2003, 302 (5649): 1364-1368.

[11] Shao QA, Gao YQ. Temperature dependence of hydrogen-bond stability in beta-hairpin structures [J]. Journal of Chemical Theory and Computation, 2010, 6(12): 3750-3760.

[12] Spassov VZ, Yan L, Flook PK. The dominant role of side-chain backbone interactions in structural realization of amino acid code. ChiRotor: a sidechain prediction algorithm based on side-chain backbone interactions [J]. Protein Science, 2007, 16(3): 494-506.

[13] Doig AJ, Sternberg MJ. Side-chain conformational entropy in protein folding [J]. Protein Science, 1995, 4: 2247-2251.

[14] Dunbrack RL, Cohen FE. Bayesian statistical analysis of protein side-chain rotamer preferences [J]. Protein Science, 1997, 6(8): 1661-1681.

[15] Wei Y, Nadler W, Hansmann UHE. Backbone and sidechain ordering in a small protein [J]. The Journal of Chemical Physics, 2008, 128: 025105.

[16] Hansmann UHE. Parallel tempering algorithm for conformational studies of biological molecules [J]. Chemical Physics Letters, 1997, 281: 140-150.

[17] Subramani A, DiMaggio PA, Floudas CA. Selecting high quality protein structures from diverse conformational ensembles [J]. Biophysical Journal, 2009, 97: 1728-1736.

[18] Anfinsen C. The formation and stabilization of protein structure [J]. Biochemical Journal, 1972, 128(4): 737-749.

[19] Mitsutake A, Sugita Y, Okamoto Y. Generalizedensemble algorithms for molecular simulations of biopolymers [J]. Biopolymers, 2001, 60: 96-123.

[20] Hansmann UHE, Berg BA, Neuhaus T. Recent results of multimagnetical simulations of the ising model [J]. International Journal of Modern Physics C, 1992, 3: 1155-1161.

[21] Wang FG, Landau DP. Efficient, multiple-range random walk algorithm to calculate the density of states [J]. Physical Review Letters, 2001, 86 (10): 2050-2053

[22] Marinari E, Parisi G. Simulated tempering: a new monte carlo scheme [J]. Europhysics Letters, 1992, 19: 451-458.

[23] Hesselbo B, Stinchcombe RB. Monte carlo simulation and global optimization without parameters [J]. Physical Review Letters, 1995, 74: 2151-2155.

[24] Hukushima K, Nemoto K. Exchange monte carlo method and application to spin glass simulations [J]. Journal of the Physical Society of Japan, 1996, 65: 1604-1608.

[25] Hansmann UHE. Parallel tempering algorithm for conformational studies of biological molecules [J]. Chemical Physics Letters, 1997, 281: 140-150.

[26] Sugita Y, Okamoto Y. Replica-exchange molecular dynamics method for protein folding [J]. Chemical Physics Letters, 1999, 314: 141-151.

[27] Hansmann UHE, Okamoto Y. New monte carlo algorithms for protein folding [J]. Current Opinion in Structural Biology, 1999, 9: 177-183.

[28] Lee J, Scheraga HA, Rackovsky S. New optimization method for conformational energy calculations on polypeptides: conformational space annealing [J]. Journal of Computational Chemistry, 1997, 18: 1222-1232.

[39] Pedersen JT, Moult J. Genetic algorithms for protein structure prediction [J]. Current Opinion in Structural Biology, l996, 6: 227-231.

[30] Androulakis IP, Maranas CD, Floudas CA. Prediction of oligopeptide conformations via deterministic global optimization [J]. Journal of Global Optimization, 1997, 11: 1-34.

[31] Samudrala R, Moult J. A graph-theoretic algorithm for comparative modeling of protein structure [J]. Journal of Molecular Biology, 1998, 279: 287-302. [32] Subramani A, Wei Y, Floudas CA. ASTRO-FOLD 2.0: an enhanced framework for protein structure prediction [J]. AIChE Journal, 2012, 58(5): 1619-1637.

[33] Carpio CAD. A parallel genetic algorithm for polypeptide three dimensional structure prediction. A transputer implementation [J]. Journal of Chemical Information and Computer Sciences, 1996, 36: 258-269.

[34] Kikugawa G, Apostolov R, Kamiya N, et al. Application of MDGRAPE-3, a special purpose board for molecular dynamics simulations, to periodic biomolecular systems [J]. Journal of Computational Chemistry, 2009, 30(1): 110-118.

[35] Shaw DE, Dror RO, Salmon JK, et al. Millisecondscale molecular dynamics simulations on Anton [C] // Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, 2009: 1-11.

[36] Grossman JP, Kuskin JS, Bank JA, et al. Hardware support for fine-grained event-driven computation in Anton 2 [C] // Proceedings of the Eighteenth International Conference on Architectural Supportfor Programming Languages and Operating Systems, 2013: 549-560.

[37] Alder BJ, Wainwright TE. Studies in molecular dynamics. I. general method [J]. The Journal of Chemical Physics, 1959, 31(2): 459-466.

[38] Rahman A. Correlations in the motion of atoms in liquid argon [J]. Physical Review, 1964, 136(2A): A405-A411.

[39] Zhao G, Perilla JR, Yufenyuy EL, et al. Mature HIV-1 capsid structure by cryo-electron microscopy and all-atom molecular dynamics [J]. Nature, 2013, 497(7451): 643-646.

[40] Shaw DE, Maragakis P, Lindorff-Larsen K, et al. Atomic-level characterization of the structural dynamics of proteins [J]. Science, 2010, 330(6002): 341-346.

[41] Arkhipov A, Shan Y, Das R, et al. Architecture and membrane interactions of the EGF receptor [J]. Cell, 2013, 152(3): 557-569.

[42] Dror RO, Green HF, Valant C, et al. Structural basis for modulation of a G-protein-coupled receptor by allosteric drugs [J]. Nature, 2013, 503(7475): 295-299.

[43] Price DJ, Brooks CL. A modified TIP3P water potential for simulation with Ewald summation [J]. The Journal of Chemical Physics, 2004, 121(20): 10096-10103.

[44] Cerutti DS, Case DA. Multi-level ewald: a hybrid multigrid/fast fourier transform approach to the electrostatic particle-mesh problem [J]. Journal of Chemical Theory and Computation, 2009, 6(2): 443-458.

[45] Walker RC, Crowley MF, Case DA. The implementation of a fast and accurate QM/MM potential method in Amber [J]. Journal of Computational Chemistry, 2008, 29(7): 1019-1031.

[46] 周亞芳, 江泓, 湯建光, 等. 蛋白質(zhì)磷酸化修飾在多聚谷氨酰胺疾病中的研究進(jìn)展 [J]. 中華醫(yī)學(xué)遺傳學(xué)雜志, 2008, 25(4): 414-417.

[47] Paulson HL, Bonini NM, Roth KA. Polyglutamine disease and neuronal cell death [J]. PNAS, 2000, 97: 12957-12958.

[48] Khare1 SD, Ding F, Gwanmesia KN, et al. Molecular origin of polyglutamine aggregation in neurodegenerative disease [J]. PLoS Computational Biology, 2005, 1(3): e30.

[49] Wanker EE. Protein aggregation and pathogenesis of huntington’s disease: mechanisms and correlstions [J]. The Journal of Biological Chemistry, 2000, 381: 937-942.

[50] Ogawa H, Nakano M, Watanabe H, et al. Molecular dynamics simulation study on the structural stabilities of polyglutamine peptides [J]. Computational Biology and Chemistry, 2008, 32(2): 102-110.

[51] Merlino A, Esposito L, Vitagliano L. Polyglutamine repeats and beta-Helix structure: molecular dynamics study [J]. Proteins: Structure, Function, and Bioinformatics, 2006, 63: 918-927.

[52] Perutz MF, Finch JT, Berriman J, et al. Amyloid fibers are water-filled nanotubes [J]. Proceedings of the National Academy of Sciences of the United States of America, 2002, 99(8): 5591-5595.

[53] Chopra M, Reddy AS, Abbott NL, et al. Folding of polyglutamine chains [J]. The Journal of Chemical Physics, 2008, 129: 135102.

[54] Barton S, Jacak R, Khare SD, et al. The length dependence of the polyQ-mediated protein aggregation [J]. The Journal of Biological Chemistry, 2007, 282(35): 25487-25492.

[55] Nelson R, Sawaya MR, Balbirnie M, et al. Structure of the cross-beta spine of amyloid-like fi brils [J]. Nature, 2005, 435: 773-778.

[56] Esposito L, Paladino A, Pedone C, et al. Insights into structure, stability, and toxicity of monomeric and aggregated polyglutamine models from molecular dynamics simulations [J]. Biophysical Journal, 2008, 94(10): 4031-4040.

[57] Marti-Renom MA, Stuart AC, Fiser A, et al. Comparative protein structure modeling of genes and genomes [J]. Annual Review of Biophysics and Biomolecular Structure, 2000, 29: 291-325.

[58] Bowie JU, Luthy R, Eisenberg D. A method to identify protein sequences that fold into a known three-dimensional structure [J]. Science, 1991, 253: 164-170.

[59] Jones DT, Taylor WR, Thornton JM. A new approach to protein fold recognition [J]. Nature, 1992, 358: 86-89.

[60] Wu S, Skolnick J, Zhang Y. Ab initio modeling of small proteins by iterative TASSER simulations [J]. BMC Biology, 2007, 5: 17.

[61] Das R, Baker D. Macromolecular modeling with rosetta [J]. Annual Review of Biochemistry, 2008, 77(1): 363-382.

[62] Kryshtafovych A, Fidelis K, Moult J. CASP10 results compared to those of previous CASP experiments [J]. Proteins, 2014, 82(Suppl 2): 164-174.

[63] Moult J, Fidelis K, Kryshtafovych A, et al. Critical assessment of methods of protein structure prediction (CASP)-round x [J]. Proteins, 2014, 82(Suppl 2): 1-6.

[64] Mirjalili V, Noyes K, Feig M. Physics based protein structure refinement through multiple molecular dynamics trajectories and structure averaging [J]. Proteins, 2014, 82(Suppl 2): 196-207.

Advances on Computer Simulations of Protein Folding

WEI Yanjie ZHANG Huiling HUANG Qingsheng

( Shenzhen Key Lab for High Performance Data Mining, Center for High Performance Computing, Shenzhen Institutes of Advanced Technology,Chinese Academy of Sciences,Shenzhen 518055,China )

Protein misfolding leads to many diseases, such as Alzhemizer disease and polyQ diseases, etc. Protein folding is vital for the mechanism study of folding-related diseases. In this paper, several aspects in protein folding were reviewed, including the side chain ordering, computational algorithms for protein folding, folding diseases, molecular dynamics simulation and structure prediction.

protein folding; parallel monte carlo method; protein structure prediction; molecular dynamics simulation; protein folding disease

TG 156

A

2013-11-29

國(guó)家自然科學(xué)基金項(xiàng)目(11204342)；深圳市科創(chuàng)委項(xiàng)(JCYJ20120615140912201)；深圳市孔雀計(jì)劃項(xiàng)目(KQCX20130628112914299)。

魏彥杰(通訊作者)，博士，研究方向?yàn)橛?jì)算生物學(xué)，E-mail：yj.wei@siat.ac.cn；張慧玲，碩士，研究方向?yàn)樯镄畔W(xué)和計(jì)算生物學(xué)；黃慶生，博士，研究方向?yàn)樯镄畔W(xué)和計(jì)算生物學(xué)。