基于HNP模型的強化學習狀態空間表示方法

吳宏杰 韓佳妍 楊 茹 陸衛忠 傅啟明*

1(蘇州科技大學電子與信息工程學院 江蘇 蘇州 215009)2(蘇州大學江蘇省計算機信息處理技術重點實驗室 江蘇 蘇州 215006)

0 引 言

蛋白質折疊構形預測問題是現如今生物信息學領域中具有高挑戰難度的問題之一，清楚蛋白質折疊構形對生物學發展和人們生活都有重大意義[1]。生物信息學假說認為自然界存在的蛋白質形態是最穩定的形態(即自由能最低)[2]。蛋白質分子的折疊過程是指從蛋白質分子從一般狀態變化到基態的復雜過程，它能使我們了解氨基酸序列是如何決定蛋白質分子結構，預測其結構及結構所表現出來的蛋白質分子的性能[3]。目前，國內外的相關研究提出了許多可計算的理論模型,如HP模型[4-6]、AB模型[7-8]和HNP模型[9-12]等。為了能夠達到簡化對序列空間搜索的目的，很多基于蛋白質設計方法應用于最簡單的氨基酸兩態分類的HP模型上并且取得了成功[5]。但是進一步研究表明，在很多情況下，簡單HP模型存在一些問題，Shakhnovich[13]指出，當蛋白質理論模型設計中考慮更多的氨基酸類型時，這種情況就會得到改善。為此本文將基于強化學習搜索蛋白質序列空間的方法應用到HNP模型上。

如今，強化學習已經在多個領域得到廣泛應用。例如車輛定位與識別、游戲自動化檢測和機器人仿真等[14]；在生物科學領域也發揮著重要作用，例如生物序列對比、基因組測序等[15]。在強化學習中，Agent可以收到環境的狀態反饋，并且與環境產生互動，通過得到的經驗來進行自主學習，根據獲得的獎賞反饋來找到整體最優解，而不會陷于局部最優。但是存在“維數災難”(即狀態空間的大小隨著特征數量的增加而發生指數級的增長)和收斂速度遲緩這兩個嚴重且普遍的問題[16]。本文選用強化學習中經典的Q-learning算法實現HNP晶格模型的最優結構預測。但Q-learning方法及其相關演變算法只在小的狀態空間范圍內比較有效，一旦訓練的狀態空間變大，則對計算機內存的要求較高，其存儲和計算相應的狀態動作對也會變得更加困難，這就是在利用強化學習算法計算蛋白質序列最優結構所遇到的“維數災難”問題。

目前對于“維數災難問題”，有研究者提出函數逼近的方法，它無須對狀態空間進行預處理，可以用較小的存儲空間來描述狀態空間到動作的映射關系，但收斂性無法得到保證，理論基礎還相當薄弱[16]。有研究者提出了分層強化學習，將復雜問題分解為多個小問題，通過解決多個小問題從而達到能夠解決原問題的目的[17]。如果利用分層強化學習思想解決蛋白質序列結構預測問題則會出現序列特征擺放重復的問題，不符合真實的實際蛋白質序列結構。因此，基于HNP晶格模型，本文利用強化學習算法計算序列的最低能量的基礎上研究強化學習的狀態空間。通過比較全狀態空間表示法與簡單狀態空間表示法，結合兩者的優勢，提出一種半狀態空間表示法，有效解決蛋白質序列特征增多導致狀態空間過大的問題，達到改善“維數災難”問題的目的。在相同的實驗環境下，比較這三種不同狀態空間表示方法的實驗結果，對比得出三種狀態空間表示方法的優缺點：1) 全狀態空間表示法在短序列的實驗上可以準確計算最低能量，但其狀態空間受計算機內存限制，無法進行長序列的能量計算；2) 簡單狀態空間表示法可以對長短序列都進行計算能量，但由于其狀態空間設置過小，在計算最低能量的準確度上大打折扣；3) 半狀態空間表示法有效解決了全狀態空間表示法和簡單狀態空間表示法的缺點，能夠對長序列的能量進行計算且比簡單狀態空間的準確度高。

1 HNP模型

1.1 概 念

很多學者不斷提出針對于蛋白質折疊結構問題的簡化模型，如HP模型[4-6]和AB模型[7-8]。而Miyazawa等[18]和Zhang等[19]采用統計方法計算了這20種氨基酸之間的相互作用能，根據相互作用能的大小，提出把20種氨基酸分為3類。根據相互作用能的大小，分別將氨基酸Cys、Met、Phe、Ile、Val、Trp、Tyr稱為疏水(H)基團，氨基酸Ala、Gly、Thr、His稱為中性基團(N)和Ser、Asn、Gln、Asp、Glu、Arg、Lys、Pro稱為親水(P)基團，即HNP模型[9]。

1.2 能量函數表示

根據氨基酸相互作用能的大小，將氨基酸分為3類：疏水氨基酸(用H表示)、中性氨基酸(用N表示)、親水氨基酸(用P表示)。通過氨基酸的劃分把蛋白質鏈轉化為以H、N、P構成的序列。

根據蛋白質分子的Hamiltonian 能量函數[20-21]：

(1)

式中：Si、Sj表示氨基酸分類H、N和P；ESi-Sj表示氨基酸與氨基酸之間的相互作用能量，EH-H=-5.59 RT，EH-N=-3.68 RT，EH-P=-3.07 RT，EN-N=-2.23 RT,EN-P=-1.78 RT,EP-P=-1.37 RT。rij指的是氨基酸與氨基酸之間的距離；l為假定的蛋白質分子的鍵長；δ(rij,l)是判斷氨基酸是否形成一個緊密接觸對的條件，如果rij=l，則δ(rij,l)=1，即形成一個緊密接觸對，如果rij≠l，則δ(rij,l)=0，即不能形成一個緊密接觸對。

2 基于HNP模型強化學習狀態空間表示方法

2.1 基于不同狀態空間的統一強化學習框架

強化學習(Reinforcement Learning, RL)是指從環境狀態到行為映射的學習，目的是使Agent從環境中獲得的積累獎賞值達到最大[22]。Q-learning算法是強化學習中一種獨立于模型的TD控制算法，它可以用來為任何給定的馬爾可夫決策過程(Markov Decision Process,MDP)找到最佳的行動選擇策略。Q-learning算法的最簡單形式被定義為：

(2)

Q-learning算法最大的優點就是無需環境模型，在學習過程中，該算法先建立一個Q值表，設其初值為0。當Agent采取動作at轉移至下一狀態st+1，并獲得立即獎賞rt+1，根據式(2)，更新Q值表，經過不斷迭代更新，Q值表會最終收斂到最優，根據Q值表選擇最優動作值函數可以選取到最優動作。

本文通過強化學習中的Q學習算法來實現HNP模型結構的預測。其流程框架如圖1所示。

圖1 基于HNP模型的強化學習框架

(1) 序列輸入與預處理。本文將氨基酸序列轉化為HNP序列的形式，在運算時，為了方便考慮，通過-1、0和1來對應簡化代替H(疏水)、N(中性)和P(親水)三類氨基酸。

(2) 基本元素。在強化學習系統中，除了Agent和環境之外，還包含有動作、狀態值函數和獎賞(回報函數)這3類基本元素。而馬爾可夫決策過程(MDP)則包含5個模型要素，分別為狀態(state)、動作(action)、策略(policy)、獎勵(reward)和回報(return)[23]。

(3) 結果輸出。通過理論可以得知，如果所有的狀態動作對能夠不斷地被更新，那么Q值表能夠被收斂到最佳結果，Agent可以根據收斂到的Q值更新表選擇最優動作，從而得出HNP模型的最優結構預測。

2.2 狀態空間表示方法

在強化學習中，狀態是指學習器發現自身當前所處的情況：對于Agent而言，狀態可以是其所有執行器的位置和角度，以及其地圖位置、電池狀態、其試圖發現的目標等[24]。在馬爾可夫系統中，所有這些變量都被賦予學習器。狀態也可以被定義為學習器可用于發現其當前狀況的所有信息的總結。這些信息包含在它所經歷動作和取樣的歷史中，意味著完整的歷史構成完美的狀態，這也是衡量一個狀態表示的標準[25]。

強化學習中的Q-learning算法是一種以馬爾可夫決策過程為理論基礎的算法。馬爾可夫性是指Agent在當前狀態s下采取當前動作a達到下一狀態st+1以及得到的立即獎賞只與當前狀態s和當前動作a有關，與之前的狀態和動作都無關。也就是說，Agent和環境在一個離散時間序列(t=0,1,2,…)的每一步中都進行交互，在每個時間步t，Agent都能得到若干環境狀態st∈S，其中S是所有可能狀態的集合。

本文將蛋白質序列中的氨基酸可能擺放的位置看作一個狀態，想要計算出蛋白質序列的最低能量也就是需要羅列出Agent采取每一步動作所有可能的狀態的集合。而當狀態空間不斷增大時，Agent與環境交互的時間越來越長，其存儲和計算相應的狀態動作對也會變得更加困難，對計算機內存要求更高。需要簡化這個歷史，而不丟失信息。因此，基于HNP蛋白質模型，對以下三種狀態空間設置方法進行研究。

(1) 全狀態空間設置方法。對于一個長度為n的蛋白質序列，在二維折疊空間中，有四個可能的折疊方向，分別為：L(左)、U(上)、R(右)和D(下)，即A=[L,U,R,D]。如圖1(a)所示，首先隨機固定第一個氨基酸的狀態，后續的每個氨基酸可能的狀態是上一個氨基酸4個可能方向L(左)、U(上)、R(右)和D(下)的集合，也就是說后一個氨基酸的所有可能狀態個數是前一個氨基酸狀態個數的4倍，狀態轉移表如表1所示。

表1 全狀態空間轉移表

續表1

隨機固定第一個氨基酸的狀態s1，第一個氨基酸狀態個數為1。s1采取四個可能的折疊方向，轉移到的第二個氨基酸的可能狀態{s2,s3,s4,s5}，則第二個氨基酸所有可能狀態個數為4。如果s2變為當前狀態時再轉移到第三個氨基酸的可能狀態{s6,s7,s8,s9}；同理，如果s3、s4、s5變為當前狀態時也轉移到第三個氨基酸的可能狀態{s10,s11,s12,s13}、{s14,s15,s16,s17}、{s18,s19,s20,s21}，則第三個氨基酸所有可能狀態個數為4×4=16。以此類推，將整個狀態集S的狀態總數視為一個初值為1、公比為4的等比數列求和：

(2) 半狀態空間設置方法。對于全狀態空間過大，簡單狀態空間略小的問題，本文提出一種半狀態空間設置方法，如圖1(b)所示，如果當前氨基酸的上一狀態采取的動作相同，則采取所有可能的動作轉移到相同的下一狀態，相當于固定了16個狀態空間，狀態轉移表如表2所示。

表2 半狀態空間轉移表

續表2

隨機固定第一個氨基酸的狀態s1，則第一個氨基酸狀態個數為1。s1采取四個可能的折疊方向，轉移到的第二個氨基酸的可能狀態{s2,s3,s4,s5}，則第二個氨基酸所有可能狀態個數為4。如果s2變為當前狀態時再轉移到第三個氨基酸的可能狀態{s6,s7,s8,s9}；同理，如果s3、s4、s5變為當前狀態時也轉移到第三個氨基酸的可能狀態{s10,s11,s12,s13}、{s14,s15,s16,s17}、{s18,s19,s20,s21}，則第三個氨基酸所有可能狀態個數為16。從第三個氨基酸開始，考慮當前動作和上一動作轉移到下一氨基酸的可能狀態。即如果s6作為當前狀態，上一動作為L，采取當前動作為L，則轉移到s22；采取當前動作為R，則轉移到s23；采取當前動作為U，則轉移到s24；采取當前動作為D，則轉移到s25。如果s7作為當前狀態，上一動作為R，采取當前動作為L，則轉移到s26；采取當前動作為R，則轉移到s27；采取當前動作為U，則轉移到s28；采取當前動作為D，則轉移到s29。如果s8作為當前狀態，上一動作為U，采取當前動作為L，則轉移到s30；采取當前動作為R，則轉移到s31；采取當前動作為U，則轉移到s32；采取當前動作為D，則轉移到s33。如果s9作為當前狀態，上一動作為D，采取當前動作為L，則轉移到s34；采取當前動作為R，則轉移到s35；采取當前動作為U，則轉移到s36；采取當前動作為D，則轉移到s37。如果s10作為當前狀態，上一動作為L，采取當前動作為L，則轉移到s22；采取當前動作為R，則轉移到s23；采取當前動作為U，則轉移到s24；采取當前動作為D，則轉移到s25。以此類推，從第三個氨基酸開始，每個氨基酸的所有可能狀態都為16，則整個狀態集S的狀態總數為：

1+41+42+42…+42=1+4+16×(n-2)=16n-27

所以一個蛋白質分子的半狀態集S={s1,s2,…,s16n-27}。

(3) 簡單狀態空間設置方法。全狀態空間的設置方法全部采取四個可能的折疊方向，全部羅列導致狀態空間過大，其存儲和計算序列的狀態動作對會比較困難。因此本文提出第一步隨機固定第一個氨基酸的狀態，從第二步開始每一個氨基酸可能的狀態都固定為四個狀態空間，圖形表示如圖1(c)所示，狀態轉移表如表3。隨機固定第一個氨基酸的狀態s1，則第一個氨基酸狀態個數為1。s1選擇四個可能的折疊方向轉移到下一氨基酸的可能狀態{s2,s3,s4,s5}，則第二個氨基酸所有可能狀態個數為4。如果s2變為當前狀態時再轉移到第三個氨基酸的所有可能狀態{s6,s7,s8,s9}；同理，如果s3、s4、s5變為當前狀態時也轉移到第三個氨基酸的所有可能狀態{s6,s7,s8,s9}，則第三個氨基酸所有可能狀態個數也為4。如果s6變為當前狀態時再轉移到下一可能狀態{s10,s11,s12,s13}；同理，如果s7、s8、s9變為當前狀態時也轉移到下一可能狀態{s10,s11,s12,s13}，則第四個氨基酸所有可能狀態個數也為4。由此類推，把整個狀態集S視為第一個氨基酸的狀態加上后面n-1個氨基酸狀態所有可能狀態的集合；狀態總數為：

1+41+41+…+41=1+4×(n-1)=4n-3

所以一個蛋白質分子的簡單狀態集S={s1,s2,s3,…,s4n-3}。

表3 簡單狀態空間轉移表

續表3

綜上所述，可以得出三種不同的狀態空間設置表示的總個數。為了便于觀察以及方便進行對比，將三個方法的狀態空間集記錄列于表4中。

表4 三種狀態集表示

2.3 獎賞函數設置與策略選擇

在強化學習中，獎賞函數把環境中感知到的狀態動作對映射為一個單獨的數字，即獎賞，表示該狀態的內在需求。Agent的唯一目標是在長期運行過程中實現獎賞最大化。本文將強化學習思想與HNP模型相結合，為了能夠找到能量最低的結構模型，也為了能夠讓Agent更好地訓練出所需要的理想化結構，本文選擇根據上述所說的Hamiltonian能量函數計算公式，將氨基酸與氨基酸之間的相互能量作為獎賞函數。因此，獎賞函數的設置如下：

(1) 當序列未完全擺放好時，其獎賞設置為0。

(2) 當序列擺放到最后一個氨基酸時，判斷氨基酸與氨基酸之間是否滿足形成緊密接觸對的條件。如果滿足，則相應的獎賞值為氨基酸之間的相互作用能的絕對值|E|；如果不滿足，則獎賞值為0。

MDP策略是按狀態給出的，動作的條件概率分布在強化學習下屬于隨機性策略。在算法的訓練階段，本文采用“探索-貪心(ε-greedy)策略”進行動作選擇。ε-greedy策略既考慮當前回報，又兼顧將來回報，算法中隨機采取概率ε(0<ε<1)選取下一個動作，利用剩下1-ε的概率去采取最優動作，從而不斷計算更新Q值表。

2.4 訓練算法

流程如算法1所示，其中：S表示狀態設置，Q表示狀態動作Q值矩陣，at表示動作，R為獎賞函數，st為當前狀態，st+1為下一狀態，terminal為終點狀態,n為當前狀態數，N表示狀態總數。

算法1HNP模型訓練算法

1. //設置全狀態空間

2. state_transfer=0

3. for n in range(np.int((4 ** N-1)/3)):

4. state_transfer+=1

5. Or //設置半狀態空間

6. state_transfer=list(range(2，22))

7. for n in range(4，N+1):

8. state_transfer+=list(range((n*16-42),n*16-26))*4

9. Or //設置簡單狀態空間

10. state_transfer=list(range(2，6))

11. for n in range(2，N):

12. state_transfer+=list(range((n*4-2),n*4+2))*4

13. Initialize:st,Q(s,a)

14. Repeat (for each episode){

15. Initializest

16. While(st!=ternimal){

17.at← Q(ε-greedy)

18. Take actionat, returnR,st+1

20.s←st+1

21. }

22. }

3 實驗與結果分析

3.1 實驗參數選擇

為了保證在實驗過程中的客觀性以及值函數的準確性和收斂性，需要設置一定參數來進行實驗。實驗中參數有學習參數γ、步長參數α和探索概率ε。學習參數γ又稱折扣率，決定著未來獎賞的當前價值。在強化學習問題中，一般取學習參數γ=0.9，使得Agent在學習探索過程中能夠找到全局最優解而不陷于局部最優解。步長參數和探索概率共同影響著值函數的準確性和收斂性。步長參數過大，會引起值函數震蕩，難以收斂；過小則使得值函數收斂速度緩慢。探索概率過大，Agent為尋求更優解會不斷更新值函數，不易收斂；過小則缺少探索信息，值函數收斂緩慢。因此需要選擇最優參數以達到收斂效果最佳的目的。

本文在設置學習參數γ=0.9的前提下，選擇不同的步長參數α和探索概率ε進行實驗。以序列HHPNHHNNHH(長度為10)作為實驗對象，進行5輪實驗，訓練迭代次數為30萬次，改變步長參數α和探索概率ε，多次實驗計算蛋白質序列在不同參數下達到最低能量所需的平均次數，從而比較得出最優參數。得出數據如表5所示。

表5 三種狀態空間設置方法在不同參數下序列

由表5所示，在探索概率ε數值相同的情況下，隨著步長參數α的增大，序列達到最低能量所需平均次數逐漸下降，當α=0.9的時候，實驗數據趨于收斂；在步長參數α數值相同的情況下，探索概率ε=0.5測試出的所需次數較ε=0.2和ε=0.8更優。綜上所述，本文選取α=0.9、ε=0.5作為最優參數。

3.2 對比實驗

根據三種狀態空間表示方法設置的不同對蛋白質序列進行模型訓練。本文在PIR數據庫中隨機選取了11條長度不同的序列作為實驗對象，將其轉化為HNP序列，將序列集信息記錄于表6。為避免偶然性，按照三種狀態空間設置的方法對每條序列進行5輪測試實驗，每輪訓練迭代次數設置為30萬次，保證其他參數條件以及實驗環境不變，記錄在迭代次數內每條序列可以計算出的最低能量并統計于表6，其中未能測出序列最低能量的結果以“-”表示。

表6 三種狀態空間方法測得HNP序列的最低能量 單位：RT

由表6可以看出，當序列長度較短時，以前面5條序列為例，基于全狀態空間方法和半狀態空間方法的實驗都可以測出序列的最低能量；而基于簡單狀態空間方法的實驗則不能測出序列3和序列5的最低能量，相比于另外兩種方法在計算序列最低能量的準確度上降低了40百分點。此外，基于簡單狀態空間方法的實驗在訓練序列2和序列4的實驗中不能每次都計算出最低能量，相比于基于全狀態空間方法的實驗在五輪訓練實驗的準確度上平均降低了12百分點；基于半狀態空間方法的實驗在訓練序列3的實驗中不能每次都計算出最低能量，相比于基于全狀態空間方法的實驗在五輪訓練實驗的準確度上平均降低了4百分點，相比于基于簡單狀態空間方法的實驗在五輪訓練實驗的準確度上平均提高了8百分點。當序列長度逐漸增大時，以后面6條序列為例，基于全狀態空間方法的實驗則會出現狀態空間過大的情況，導致在相同實驗環境下無法計算后面6條長序列的最低能量；基于半狀態空間方法的實驗相比于基于簡單狀態空間方法的實驗能計算出更低的能量，尤其對于序列6、7、9、10和11，基于半狀態空間方法的實驗比基于簡單狀態空間方法的實驗在每輪訓練實驗中都計算出更低的能量。由此可見，對于長序列而言，半狀態空間方法比全狀態方法和簡單狀態空間方法能找到更符合實際蛋白質序列的最低能量。

綜上所述，半狀態空間方法和簡單狀態空間方法有效解決全狀態空間無法計算長序列的缺點。簡單狀態空間較全狀態空間有3條序列預測出更低能量，半狀態空間較全狀態空間方法全部6條長序列都預測出更低能量，且半狀態空間預測的能量平均值較簡單狀態空間降低了9.83百分點。

4 結 語

本文以強化學習算法為核心，對蛋白質HNP模型結構進行訓練模型，計算最接近真實的蛋白質序列結構。針對強化學習中的狀態空間問題，對比研究了三種不同的狀態空間設置方法。全狀態空間表示方法全面地列出了序列所有的可能狀態空間，能夠較為準確地計算出最接近真實的蛋白質序列結構，但每當序列長度增加時，狀態空間也會以指數級增長，導致狀態空間過大，其存儲和計算狀態動作對較為困難，對計算機內存的要求也變高。半狀態空間表示方法和簡單狀態空間表示方法在全狀態空間表示方法的基礎上將其狀態空間進行簡化，縮小了狀態空間，可以計算長度更長的序列，但也降低了一定的準確度。此外，半狀態空間表示方法相較于簡單狀態空間表示方法在計算蛋白質序列最低能量的準確率上更具優勢。

本文是強化學習在生物信息領域新的嘗試。目前還存在需要解決的問題，如無法在確保準確度的前提下合理縮小氨基酸序列的狀態空間。所以，在研究HNP模型的蛋白質分子構象和折疊問題中，需要提出更進一步的算法改進。