禁忌搜索算法在蛋白質結構預測中應用

摘要：一種氨基酸序列只可能有一種蛋白質結構，所以在蛋白質理論預測中，正確定義能量函數、精確選用的計算機搜尋算法來尋找能量最低值，是蛋白質結構預測的關鍵。基于此，本文以兩兩殘基之間距離分布和二面角分布符合玻爾茲曼定理，提出了一種抽象的蛋白質三維結構連續物理數學模型。然后應用了禁忌搜索算法很好的計算了牛胰島素B(D)主鏈走向；比較計算了氨基酸序列最低能量的全局最優點。

關鍵詞：禁忌搜索算法；蛋白質結構；預測；距離分布；二面角分布

中圖分類號:TP311.13文獻標識碼:A文章編號：1009-3044(2008)27-2101-03

The Application of Taboo Search Algorithm in Protein Structure Prediction

WANG Jian

(Department of Computer and Information Science，Neijiang Teachers College， Neijiang 641112， China)

Abstract: An amino acid alignment has only one protein structure， so in the theory of protein prediction， it is crucial to define energy function and to find minimum energy by using computer search algorithm. Wherein， this paper proposes the corresponding constant physical mathematical model of the 3D space protein according to the distance distribution and dihedral angel distribution conforming with Boltzmann distribution. Furthermore， with the application of taboo search algorithm， stretches of the B (D) main chain of the Bovine Insulin-like has been calculated and the most proper point of the minimum energy in amino acid alignment has been calculated.

Key words:taboo search algorithm;protein structure;prediction;distance distribution;dihedral angel distribution

1 引言

蛋白質三級結構是體現其生化功能和細胞功能的基礎[1]。蛋白質預測的目標是三級結構和功能預測，最終實現分子設計和藥物合成。蛋白質的空間結構極其復雜，現在有兩種方法確定：一種是實驗測量，包括用X射線衍射和核磁共振成像；另一種是理論預測，利用計算機根據理論和已知的氨基酸序列等信息來預測，方法包括同源結構模擬、折疊辨識模擬和基于第一性原理的從頭計算[3-4]。對于蛋白質結構的實驗測定十分費時費力，如今隨著技術的進步，實驗測蛋白質結構的時間和花費已經大大地減少了，但測定一個蛋白質結構的平均費用也在100萬美元左右。而對于計算機的預測，雖然對同源蛋白結構預測得到了很好的結果，但對于無同源信息的蛋白質預測是非常困難的，無同源信息的蛋白質預測難以從預測模型中提取領域相關知識，人們對蛋白質的構象機理認識有限，以至于蛋白質結構預測準確率有待提高。由于一種氨基酸序列只可能有一種蛋白質結構，因此建立一種好的蛋白質能量模型和選擇一種好的學習算法是提高蛋白質預測精確率關鍵所在[2]。針對這些困難與挑戰，分析現有預測方法基礎上，從蛋白質氨基酸殘基間距離分布和二面角分布符合玻爾茲曼定理，研究提出了蛋白質距離和二面角物理數學模型，并將禁忌搜索算法[5-7]應用于蛋白質結構預測問題。

2 物理數學模型建立

從氨基酸序列預測三級結構是蛋白質結構預測的最終目標。建立預測模型和模型的求解是很關鍵的，預測模型的核心是目標函數，它至少應該能夠區分蛋白質天然構象和其它構象，通常根據目標函數建立的基礎分為基于物理理論的經驗勢能函數和基于統計方法的平均勢能函數。本文根據蛋白質分子在溶液中微觀狀態的分布符合玻爾茲曼定理，所以把兩兩氨基酸之間的距離分布和二面角分布用玻爾茲曼定理來描述。

2.1 氨基酸殘基之間的距離分布

玻爾茲曼定理描述了一個系統的自由能由構成它的各個粒子的內能及這些粒子在不同微觀狀態的分布形成的熵所決定[5]。粒子的分布不僅按速率區間分布，還應按位置區間分布，考慮到蛋白質分子各殘基在自由能最低狀態下形成一個穩定結構，速度的分布不是主要的，所以從粗粒度上簡化兩兩氨基酸殘基之間的作用，以兩兩氨基酸殘基之間的距離來設計蛋白質三級結構物理數學模型。在系統的N個粒子中，分布在某能級i上的粒子數正比于該能級的多重度gi與玻耳茲曼因子exp[-Ei/kT]的乘積，寫成數學表達式：

（1）

因為離子分配函數是粒子在各能級上有效容量的總和，引入分配函數q。，（2）式改寫為

（2）

（2）式中q=N/λ，為了表示兩個殘基ai和aj的一維序列間距φ(φ=i-j)的影響，把各能量級考慮為與相鄰氨基酸殘基的距離r所影響的，所以（2）式可改寫為

（3）

在（3）式中αi，αj分別表示20種氨基酸之中一種，i，j表示兩個氨基酸殘基的序列次序， pφαi，αj(r)表示當αi，αj之間距離為r時的概率密度，qφαi，αj表示分配函數，Eφαi，αj(r)表示當概率密度為qφαi，αj時的能量大小，k表示玻耳茲曼常數，T表示絕對溫度。為便于計算絕對能量Eφαi，αj(r)的值，在已知一確定序列情況下分配函數qφαi，αj表示蛋白質構象也確定。所以(4)式用能量函數表示為：

（4）

2.2 二面角分布

在描述肽鍵的剛性和可變程度時，運用了兩面角定量地描述α碳原子和肽平面間單鍵的旋轉，運用α碳原子和氨基一側肽平面間兩面角Ф和α碳原子和羧基一側的肽平面間兩面角Ψ如圖1所示[7]。為便于研究假設氨基酸二面角形成的能量值與兩兩氨基酸之間的距離分布相同，建立一個類似的二面角能量函數

（5）

2.3 組合能量函數

在蛋白質結構中，主鏈二面角（Ф，Ψ）分布符合拉氏構象圖，兩兩氨基酸之間的距離分布也有規律可循。因此，可利用（4）、（5）式組合來構建基于距離和二面角的能量函數。

（6）

3 禁忌搜索算法設計

3.1 禁忌搜索算法的原理和實現步驟

利用禁忌搜索算法求解組合優化問題時，首先按照隨機方法產生一個初始解作為當前解，然后在當前解的鄰域中搜索若干個解，取其中的最好解作為新的當前解。為了避免對已搜索過的局部最優解的重復，禁忌搜索算法使用禁忌表記錄已搜索的局部最優解的歷史信息，這使得算法可在一定程度上避開局部最優點，從而開辟新的搜索區域。

用禁忌搜索算法求解組合優化問題時，其實現步驟為（以目標函數求最小為例）：

第一步：選定一個初始解xnow；令禁忌表H=φ；

第二步：若滿足終止準則，轉第四步；否則，在xnow的鄰域N（xnow）中選出滿足禁忌要求的候選集C-N（xnow），轉第三步；

第三步：在C-N（xnow）中選一個評價值最好的解xbest，令xnow = xbest，更新禁忌表H，轉第二步；

第四步：輸出計算結果，停止。

3.2 禁忌搜索算法設計

基于上述蛋白質模型抽象表示方法，結合蛋白質計算機求解問題的特點，作者構造了求解該問題的一種新的禁忌搜索算法，其算法策略如下。

1）解的評價方法。用禁忌搜索算法求解組合優化問題時，需要對解進行評價，使算法在迭代過程中，不斷搜索到質量更優的解。根據蛋白質模型的描述，對于某個解要判定其優劣，首先要得到距離分布有解這一方案，然后要判斷二面角分布是否滿足問題的約束條件，同時計算距離分布和二面角分布的目標函數值，在滿足問題的約束條件的前提下，其目標函數值越優，則解的質量越高。采用距離分布表示方法產生的解所確定的模型，隱含能夠滿足二面角分布這一約束條件。對于某個解，設其對應的距離分布與二面角分配之差為M，其目標函數值為Z，可將M看成該解對應的無效長程數，設對每條無效長程的懲罰權重為Pw（該權重可根據目標函數的取值范圍取一個相對較大的正數），則該解的評價值E可用公式（7）計算。

2）鄰域操作方法。禁忌搜索算法是一種基于鄰域搜索技術的算法，確定鄰域操作方法是構造該算法的一個重要步驟。本文采用兩交換法實施鄰域操作。該方法是指隨機選擇解中的兩個元素，并交換其值的鄰域操作方法。

3）禁忌對象的確定。禁忌對象是指禁忌表中被禁的那些局部最優解。本文將每次迭代得到的最好解作為禁忌對象放入禁忌表中。

4）禁忌長度的確定。禁忌長度是指被禁對象不允許被選取的迭代步數。本文取禁忌長度為一個常數，其值應根據問題的規模確定。

5）候選集合的確定。本文將從當前解的鄰域中隨機選擇若干個鄰居作為候選集合。

6）終止準則的確定。本文采用迭代指定步數的終止準則。

4 實驗結果及分析

實驗1：本文選用PDB數據庫中牛胰島素B(D)鏈作為計算對象。牛胰島素含有4條鏈，空間結構呈典型的全α螺旋結構域，其中B(D)鏈為為一個α螺旋結構域，由25個氨基酸殘基組成，其氨基酸序列為：H-HPE-VαL- ASN-GLN-HIS-LEU-CYS-GLY-SER-HIS-LEU-VAL-GLU-ALA-LEU-TYR-LEU-AL-CYS-GLY-GLU-ARG-GLY- PHE- PHE-OH。從第8個氨基酸殘基GLY到第22個氨基酸殘基ARG為一個完整α螺旋結構。對能量變化起作用的變量有91個，能量表面的局部極小點超過390個。利用禁忌搜索算法，7次運算3次收斂到同一極小點，極小值為-133.421kal/mol。計算結果表明B(D)主鏈走向與其空間結構呈典型的全α螺旋結構域這一結論相符合。

實驗2：選用斐波納挈序列作為計算對象，設置如下：S0=A，S1=B，Si+1=Si-1Si，其中為連接算子。Si由斐波納挈序列給出，那么在序列中疏水性殘基A被孤立，親水性殘基B成對或單個被孤立。通過與PERM方法最優能量值Eperm及共軛梯度法優能量值Emin比較[9，10]，說明模型設計和算法選取是可行的，比較結果數據見表1。

5 結論

從上面的計算結果上看，本文從距離和二面角出發，設計的物理數學模型能夠表征蛋白質的特征，反映蛋白質結構的一些性質。禁忌搜索算法能較很好計算最低能量的全局最優點。

本文作者創新點以兩兩殘基之間距離分布和二面角分布符合玻爾茲曼定理，提出了一種抽象的蛋白質三維結構連續物理數學模型。然后應用了禁忌搜索算法很好的計算了牛胰島素B(D)主鏈走向；比較計算了氨基酸序列最低能量的全局最優點。總的來說：一是要正確定義能量函數，二是選用了一種精確的計算機搜尋引擎來尋找能量最低值。

參考文獻：

[1] Afinsen C B.Principles that govern the folding of protein chains[J].Science，1973，181:223-230.

[2] 王勇獻，蛋白質二級結構預測的模型和方法研究，博士學位論文，國防科技大學，2004.

[3] Jin LX and Tang HW.An improvement of simulated annealing algorithm and its application to protein structure prediction. Systems Engineering Theory and Pracitice，2002，22(9):92-96.

[4] Ni HC，Wang YF and Shi DH.Application of genetic algorithms to protein in structure prediction.Journal of Shanghai University (Natural Science)，2001，7(3):225-230.

[5] Glover F.Tabu search:part I.Orsa.Journal on Computing，1989，1:190-206.

[6] Glover F.Tabu search:part II.Orsa.Journal on Computing， 1990，2:4-32.

[7] Hertz A and de Werra D.The tabu search metaheuristic:how we use it.Annals Mathematics and Artifical Intelligence，1990，I:111-121.

[8] 王光信，劉澄凡，張積樹，物理化學（第二版）[M].北京：業出版社.2001:188-207.

[9] Stillinger FH.Collective aspects of protein folding illustrated by a toy model.Physical Review E， 1995， 52(32):2872-2877.

[10] Hsiao-Ping Hsu，Vishal Mehra and Peter Grassberger.Structure optimization in an off-lattice protein model.Physical Review E，2003，68(3):037703.

注：“本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文。”